diff --git a/src/arome/aux/mode_argslist_ll.F90 b/src/arome/aux/mode_argslist_ll.F90
index 42eaf543e67b98ba9d5b03c4ac7c5b6f84aaf79c..ed5c4f1796c8fdb592c31a474a2b92c6b46d1b71 100644
--- a/src/arome/aux/mode_argslist_ll.F90
+++ b/src/arome/aux/mode_argslist_ll.F90
@@ -5,8 +5,8 @@ CONTAINS
 !
   SUBROUTINE CLEANLIST_ll(TPLIST)
 IMPLICIT NONE
-
     TYPE(LIST_ll),  POINTER :: TPLIST ! List of fields
+    CALL ABORT
   END SUBROUTINE CLEANLIST_ll
 !
   SUBROUTINE ADD3DFIELD_ll(TPLIST, PFIELD, HNAME)
@@ -17,5 +17,6 @@ IMPLICIT NONE
   !                                              of fields
     character(len=*), intent(in) :: HNAME ! Name of the field to be added
   !
-  END SUBROUTINE ADD3DFIELD_ll
+   CALL ABORT  
+END SUBROUTINE ADD3DFIELD_ll
 END MODULE MODE_ARGSLIST_ll
diff --git a/src/arome/aux/mode_io_field_write.F90 b/src/arome/aux/mode_io_field_write.F90
index 39e6f7fec9a7603ad617ff3921ddc97eb788ae5a..e1ea4aef4651b7a11ce0bd10319e08bfb59c1edb 100644
--- a/src/arome/aux/mode_io_field_write.F90
+++ b/src/arome/aux/mode_io_field_write.F90
@@ -10,6 +10,7 @@ SUBROUTINE IO_FIELD_WRITE(TPFILE,TZFIELD,PFIELD)
     TYPE(TFIELDDATA), INTENT(IN)          :: TZFIELD
     REAL, DIMENSION(:,:,:),    INTENT(IN) :: PFIELD   ! array containing the data field
     !
+    CALL ABORT
 END SUBROUTINE IO_FIELD_WRITE
 END MODULE MODE_IO_FIELD_WRITE
 
diff --git a/src/arome/aux/mode_ll.F90 b/src/arome/aux/mode_ll.F90
index d3626160c78620753d32020d3f40b40a43bc238c..a71daecfab08c9a546172b362a32210dec4b498c 100644
--- a/src/arome/aux/mode_ll.F90
+++ b/src/arome/aux/mode_ll.F90
@@ -12,12 +12,14 @@ CONTAINS
   KYOR=1+JPHEXT
   KXEND=KSIZE1-JPHEXT
   KYEND=KSIZE2-JPHEXT
+  CALL ABORT
   END SUBROUTINE GET_INDICE_ll
 
   SUBROUTINE UPDATE_HALO_ll(TPLIST, KINFO)
   USE MODD_ARGSLIST_ll, ONLY : LIST_ll  
   TYPE(LIST_ll), POINTER :: TPLIST ! pointer to the list of fields to be updated
   INTEGER                :: KINFO  ! return status
+  CALL ABORT
   END SUBROUTINE UPDATE_HALO_ll
 
 LOGICAL FUNCTION LNORTH_ll()
diff --git a/src/arome/aux/mode_sources_neg_correct.F90 b/src/arome/aux/mode_sources_neg_correct.F90
new file mode 100644
index 0000000000000000000000000000000000000000..1b49a6e2b22a64880bea7e75e5c32001a106ecda
--- /dev/null
+++ b/src/arome/aux/mode_sources_neg_correct.F90
@@ -0,0 +1,19 @@
+MODULE MODE_SOURCES_NEG_CORRECT
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE SOURCES_NEG_CORRECT(HCLOUD, HBUDNAME, KRR, PTSTEP, PPABST, &
+                              &PTHT, PRT, PRTHS, PRRS, PRSVS, PRHODJ)
+IMPLICIT NONE
+CHARACTER(LEN=*),            INTENT(IN)           :: HCLOUD   ! Kind of cloud parameterization
+CHARACTER(LEN=*),            INTENT(IN)           :: HBUDNAME ! Budget name
+INTEGER,                     INTENT(IN)           :: KRR      ! Number of moist variables
+REAL,                        INTENT(IN)           :: PTSTEP   ! Timestep
+REAL, DIMENSION(:, :, :),    INTENT(IN)           :: PPABST   ! Absolute pressure at time t
+REAL, DIMENSION(:, :, :),    INTENT(IN)           :: PTHT     ! Theta at time t
+REAL, DIMENSION(:, :, :, :), INTENT(IN)           :: PRT      ! Moist variables at time t
+REAL, DIMENSION(:, :, :),    INTENT(INOUT)        :: PRTHS    ! Source terms
+REAL, DIMENSION(:, :, :, :), INTENT(INOUT)        :: PRRS     ! Source terms
+REAL, DIMENSION(:, :, :, :), INTENT(INOUT)        :: PRSVS    ! Source terms
+REAL, DIMENSION(:, :, :),    INTENT(IN), OPTIONAL :: PRHODJ   ! Dry density * jacobian
+END SUBROUTINE SOURCES_NEG_CORRECT
+END MODULE MODE_SOURCES_NEG_CORRECT
diff --git a/src/arome/aux/shuman.F90 b/src/arome/aux/shuman.F90
index 47109f11e7a325caf881364330a5289757ce29c5..4aaa979958fc557c5501d0d06f1748012c23b79a 100644
--- a/src/arome/aux/shuman.F90
+++ b/src/arome/aux/shuman.F90
@@ -88,6 +88,7 @@ PMXF=PA
 !END DO
 !
 !PMXF(IIU,:,:)    = PMXF(2*JPHEXT,:,:)
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
 !
 !-------------------------------------------------------------------------------
 !
@@ -182,6 +183,8 @@ PMXM=PA
 !
 !PMXM(1,:,:)    = PMXM(IIU-2*JPHEXT+1,:,:)
 !
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+
 !-------------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('MXM',1,ZHOOK_HANDLE)
@@ -370,6 +373,7 @@ PMYM=PA
 !
 !PMYM(:,1,:)    = PMYM(:,IJU-2*JPHEXT+1,:)
 !
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
 !-------------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('MYM',1,ZHOOK_HANDLE)
@@ -622,6 +626,7 @@ DO JI=1,IIU-1
 END DO
 !
 !PDXF(IIU,:,:)    = PDXF(2*JPHEXT,:,:)
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
 !
 !-------------------------------------------------------------------------------
 !
@@ -713,6 +718,7 @@ END DO
 !
 PDXM(1,:,:)    =  PDXM(IIU-2*JPHEXT+1,:,:)
 !
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
 !-------------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('DXM',1,ZHOOK_HANDLE)
@@ -803,6 +809,8 @@ END DO
 !
 !PDYF(:,IJU,:)    = PDYF(:,2*JPHEXT,:)
 !
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+
 !-------------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('DYF',1,ZHOOK_HANDLE)
@@ -892,6 +900,7 @@ DO JJ=2,IJU
 END DO
 !
 PDYM(:,1,:)    =  PDYM(:,IJU-2*JPHEXT+1,:)
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/arome/ext/aro_turb_mnh.F90 b/src/arome/ext/aro_turb_mnh.F90
index 09944b25d1d3450aff54cc1a06aed263cf96cb5c..d9609464f953c8421f9a75c5ac136711125cbf2a 100644
--- a/src/arome/ext/aro_turb_mnh.F90
+++ b/src/arome/ext/aro_turb_mnh.F90
@@ -71,6 +71,7 @@
 USE MODD_CONF
 USE MODD_CST
 USE MODD_PARAMETERS
+USE MODD_IO, ONLY: TFILEDATA
 USE MODD_BUDGET, ONLY: NBUDGET_RI, TBUDGETDATA
 !
 USE MODI_TURB
@@ -175,7 +176,8 @@ TYPE(TLDDH),   INTENT(IN), TARGET      :: YDLDDH
 TYPE(TMDDH),   INTENT(IN), TARGET      :: YDMDDH
 !
 !
- TYPE(TBUDGETDATA), DIMENSION(NBUDGET_RI) :: YLBUDGET !NBUDGET_RI is the one with the highest number needed for turb
+TYPE(TBUDGETDATA), DIMENSION(NBUDGET_RI) :: YLBUDGET !NBUDGET_RI is the one with the highest number needed for turb
+TYPE(TFILEDATA) :: ZTFILE !I/O for MesoNH
 !*       0.2   Declarations of local variables :
 !
 INTEGER :: JRR,JSV       ! Loop index for the moist and scalar variables
@@ -199,8 +201,6 @@ CHARACTER(LEN=4),DIMENSION(2)  :: HLBCX, HLBCY  ! X- and Y-direc LBC
 
 INTEGER       :: ISPLIT        ! number of time-splitting
 
-LOGICAL       ::  OCLOSE_OUT   ! Conditional closure of
-                                                   ! the OUTPUT FM-file
 LOGICAL       ::  OTURB_FLX    ! switch to write the
                                ! turbulent fluxes in the syncronous FM-file
 LOGICAL       ::  OTURB_DIAG   ! switch to write some
@@ -212,11 +212,6 @@ CHARACTER(LEN=4)   ::  HTURBDIM     ! dimensionality of the
 CHARACTER(LEN=4)   ::  HTURBLEN     ! kind of mixing length
 
 REAL          ::  ZIMPL        ! degree of implicitness
-
-CHARACTER(LEN=4)   ::  HFMFILE      ! Name of the output
-                                    ! FM-file
-CHARACTER(LEN=4)   ::  HLUOUT       ! Output-listing name for
-                                    ! model n
 !
 REAL, DIMENSION(KLON,1,KLEV+2)   :: ZDXX,ZDYY,ZDZZ,ZDZX,ZDZY
                                         ! metric coefficients
@@ -265,6 +260,9 @@ IKE=KKU-JPVEXT_TURB*KKL
 ! Numero du modele si grid nestind, toujours egal a 1
 IMI=1
 
+! Fichier I/O pour MesoNH (non-utilise dans AROME)
+ZTFILE%LOPENED=.FALSE.
+
 ! Type de condition � la limite. En 1D, sans importance. A etudier en 3D.
 HLBCX(:)='CYCL'
 HLBCY(:)='CYCL'
@@ -273,9 +271,6 @@ HLBCY(:)='CYCL'
 ISPLIT=1
 
 ! pour ecriture et diagnostic dans mesoNH, � priori les switches toujours � .F.
-OCLOSE_OUT=.FALSE.
-HFMFILE=' '
-HLUOUT= ' '
 OTURB_FLX=.FALSE.
 OTURB_DIAG=.FALSE.
 
@@ -422,18 +417,18 @@ ENDDO
 
 CL=HINST_SFU
 CALL TURB (KLEV+2,1,KKL,IMI, KRR, KRRL, KRRI, HLBCX, HLBCY, ISPLIT,IMI, &
-   & OCLOSE_OUT,OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
+   & OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
    & HTURBDIM,HTURBLEN,'NONE','NONE', CL,           &
    & HMF_UPDRAFT,ZIMPL,                                    &
-   & 2*PTSTEP, 2*PTSTEP, 2*PTSTEP,                         &
-   & HFMFILE,HLUOUT,ZDXX,ZDYY,ZDZZ,ZDZX,ZDZY,ZZZ,          &
+   & 2*PTSTEP,ZTFILE,                                      &
+   & ZDXX,ZDYY,ZDZZ,ZDZX,ZDZY,ZZZ,          &
    & ZDIRCOSXW,ZDIRCOSYW,ZDIRCOSZW,ZCOSSLOPE,ZSINSLOPE,    &
    & PRHODJ,PTHVREF,PRHODREF,                              &
    & PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
    & PPABSM,PUM,PVM,PWM,PTKEM,ZSVM,PSRCM,                  &
    & PLENGTHM,PLENGTHH,MFMOIST,                            &
    & ZBL_DEPTH,ZSBL_DEPTH,                                 &
-   & PUM,PVM,PWM,ZCEI,ZCEI_MIN,ZCEI_MAX,ZCOEF_AMPL_SAT,    &
+   & ZCEI,ZCEI_MIN,ZCEI_MAX,ZCOEF_AMPL_SAT,    &
    & PTHM,ZRM, &
    & PRUS,PRVS,PRWS,PRTHS,ZRRS,ZRSVS,PRTKES_OUT,         &
    & ZHGRAD,PSIGS,                                         &
diff --git a/src/arome/gmkpack_ignored_files b/src/arome/gmkpack_ignored_files
index 958afa751bc4a6c33fbc1e340431dadccb12c60e..31a2f877240fec698d5f7271198f3b2b4fc29fd2 100644
--- a/src/arome/gmkpack_ignored_files
+++ b/src/arome/gmkpack_ignored_files
@@ -116,3 +116,14 @@ phyex/turb/tridiag_wind.F90
 phyex/turb/tridiag_thermo.F90
 phyex/turb/tridiag_tke.F90
 phyex/turb/tridiag_massflux.F90
+phyex/turb/bl89.F90
+phyex/turb/etheta.F90
+phyex/turb/emoist.F90
+phyex/turb/rmc01.F90
+phyex/turb/sbl_depth.F90
+phyex/turb/th_r_from_thl_rt_1d.F90
+phyex/turb/th_r_from_thl_rt_2d.F90
+phyex/turb/th_r_from_thl_rt_3d.F90
+phyex/turb/thl_rt_from_th_r_mf.F90
+phyex/turb/bl_depth_diag_3d.F90
+phyex/turb/bl_depth_diag_1d.F90
diff --git a/src/arome/micro/modd_cst.F90 b/src/arome/micro/modd_cst.F90
index 1f5d39b521df169340f648dcd59e73afbe32b217..86361573695350b804440998b3b87df7b7b671bf 100644
--- a/src/arome/micro/modd_cst.F90
+++ b/src/arome/micro/modd_cst.F90
@@ -96,5 +96,6 @@ INTEGER, SAVE :: NDAYSEC        ! Number of seconds in a day
 REAL,SAVE :: RDSRV              !  XRD/XRV
 REAL,SAVE :: RDSCPD             !  XRD/XCPD
 REAL,SAVE :: RINVXP00           !  1./XP00
-
+!
+REAL,SAVE     :: XMNH_EPSILON       ! minimum space with 1.0
 END MODULE MODD_CST
diff --git a/src/arome/turb/bl_depth_diag_1d.F90 b/src/arome/turb/bl_depth_diag_1d.F90
deleted file mode 100644
index c477ca8ef4893d2c19ab0f4d68b580b95de1a775..0000000000000000000000000000000000000000
--- a/src/arome/turb/bl_depth_diag_1d.F90
+++ /dev/null
@@ -1,38 +0,0 @@
-!     ######spl
-FUNCTION BL_DEPTH_DIAG_1D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-!
-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
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',0,ZHOOK_HANDLE)
-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)
-!
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',1,ZHOOK_HANDLE)
-END FUNCTION BL_DEPTH_DIAG_1D
diff --git a/src/arome/turb/compute_entr_detr.F90 b/src/arome/turb/compute_entr_detr.F90
index 16d23c24e0da948d2781d9661395ec9b4b36065e..8211b6bd81f58275b3cd8121a3b7e614f8825538 100644
--- a/src/arome/turb/compute_entr_detr.F90
+++ b/src/arome/turb/compute_entr_detr.F90
@@ -66,7 +66,7 @@ USE MODD_CST
 !
 USE MODD_CMFSHALL
 !
-USE MODI_TH_R_FROM_THL_RT_1D 
+USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D 
 
 USE MODE_THERMO
 
diff --git a/src/arome/turb/compute_updraft.F90 b/src/arome/turb/compute_updraft.F90
index b0047fc1ec66cfacbeb6093e8910ee1f81381b69..d2ea158e8dc4333cfc67c6ede04e102221342bd5 100644
--- a/src/arome/turb/compute_updraft.F90
+++ b/src/arome/turb/compute_updraft.F90
@@ -60,7 +60,7 @@ USE MODD_CMFSHALL
 USE MODD_TURB_n, ONLY : CTURBLEN
 
 USE MODI_COMPUTE_ENTR_DETR
-USE MODI_TH_R_FROM_THL_RT_1D
+USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
 USE MODI_SHUMAN_MF, ONLY: MZM_MF, MZF_MF, GZ_M_W_MF
 
 USE MODI_COMPUTE_BL89_ML
diff --git a/src/arome/turb/compute_updraft_raha.F90 b/src/arome/turb/compute_updraft_raha.F90
index 337696338c23e4075da3b4416e6b46a74bd16957..ae0ecb2e84a66868aab21337774b032372c0bf5d 100644
--- a/src/arome/turb/compute_updraft_raha.F90
+++ b/src/arome/turb/compute_updraft_raha.F90
@@ -54,7 +54,7 @@
 USE MODD_CST
 USE MODD_CMFSHALL
 
-USE MODI_TH_R_FROM_THL_RT_1D
+USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
 USE MODI_SHUMAN_MF, ONLY: MZM_MF
 
 IMPLICIT NONE
diff --git a/src/arome/turb/compute_updraft_rhcj10.F90 b/src/arome/turb/compute_updraft_rhcj10.F90
index c6d108cf8fe3d7b6ffc3a1d82615587028a2beb7..0491bebffb06a6bb71533fd40812ab9ef0b61076 100644
--- a/src/arome/turb/compute_updraft_rhcj10.F90
+++ b/src/arome/turb/compute_updraft_rhcj10.F90
@@ -53,7 +53,7 @@
 USE MODD_CST
 USE MODD_CMFSHALL
 USE MODD_TURB_n, ONLY : CTURBLEN
-USE MODI_TH_R_FROM_THL_RT_1D
+USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
 USE MODI_SHUMAN_MF, ONLY: MZF_MF, MZM_MF, GZ_M_W_MF
 
 USE MODI_COMPUTE_BL89_ML
diff --git a/src/arome/turb/modi_bl89.F90 b/src/arome/turb/modi_bl89.F90
deleted file mode 100644
index e451772994675a8640f6cce983666f77da09a8aa..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_bl89.F90
+++ /dev/null
@@ -1,22 +0,0 @@
-!     ######spl
-      MODULE MODI_BL89
-!     ################
-INTERFACE
-      SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
-!
-INTEGER,                  INTENT(IN)  :: KKA 
-INTEGER,                  INTENT(IN)  :: KKU 
-INTEGER,                  INTENT(IN)  :: KKL 
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PZZ
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PDZZ
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTHVREF
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTHLM
-INTEGER,                  INTENT(IN)  :: KRR
-REAL, DIMENSION(:,:,:,:), INTENT(IN)  :: PRM
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTKEM
-REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PSHEAR
-REAL, DIMENSION(:,:,:),   INTENT(OUT) :: PLM
-
-END SUBROUTINE BL89
-END INTERFACE
-END MODULE MODI_BL89
diff --git a/src/arome/turb/modi_bl_depth_diag.F90 b/src/arome/turb/modi_bl_depth_diag.F90
deleted file mode 100644
index 778761f9c327fd3efa4447d3c465558cabcf8436..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_bl_depth_diag.F90
+++ /dev/null
@@ -1,36 +0,0 @@
-!     ######spl
-     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
diff --git a/src/arome/turb/modi_bl_depth_diag_3d.F90 b/src/arome/turb/modi_bl_depth_diag_3d.F90
deleted file mode 100644
index d2ea0c80dc34a1a8386177160e2374a1f02670c9..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_bl_depth_diag_3d.F90
+++ /dev/null
@@ -1,24 +0,0 @@
-!     ######spl
-     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
diff --git a/src/arome/turb/modi_emoist.F90 b/src/arome/turb/modi_emoist.F90
deleted file mode 100644
index 57e762b71b4672fa94ececaa77b7ee7e028333ed..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_emoist.F90
+++ /dev/null
@@ -1,26 +0,0 @@
-!     ######spl
-MODULE MODI_EMOIST
-!#################
-!
-INTERFACE
-!
-FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
-!
-INTEGER                              :: KRR        ! number of moist var.
-INTEGER                              :: KRRI       ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM    ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM      ! Mixing ratios, where
-!                                    PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PAMOIST   ! Amoist
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PEMOIST ! result
-!
-END FUNCTION EMOIST
-!
-END INTERFACE
-!
-END MODULE MODI_EMOIST
diff --git a/src/arome/turb/modi_etheta.F90 b/src/arome/turb/modi_etheta.F90
deleted file mode 100644
index bc08b8a8944008eb082617954334cd06d697c10e..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_etheta.F90
+++ /dev/null
@@ -1,28 +0,0 @@
-!     ######spl
-MODULE MODI_ETHETA
-!#################
-!
-INTERFACE
-!
-FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
-!
-INTEGER                              :: KRR          ! number of moist var.
-INTEGER                              :: KRRI         ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM     ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM       ! Mixing ratios, where
-!                                      PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PATHETA   ! Atheta
-!                                                    
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PETHETA ! result
-!
-!
-END FUNCTION ETHETA
-!
-END INTERFACE
-!
-END MODULE MODI_ETHETA
diff --git a/src/arome/turb/modi_prandtl.F90 b/src/arome/turb/modi_prandtl.F90
deleted file mode 100644
index 74fb802bfb827554494503afc12486fc191ca5d3..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_prandtl.F90
+++ /dev/null
@@ -1,72 +0,0 @@
-!     ######spl
-     MODULE MODI_PRANDTL
-!    ###################
-!
-INTERFACE
-!
-      SUBROUTINE PRANDTL(KKA,KKU,KKL,KRR,KRRI,OCLOSE_OUT,OTURB_DIAG,&
-                         HTURBDIM,                             &
-                         HFMFILE,HLUOUT,                       &
-                         PDXX,PDYY,PDZZ,PDZX,PDZY,             &
-                         PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
-                         PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
-                         PREDTH1,PREDR1,                       &
-                         PRED2TH3, PRED2R3, PRED2THR3,         &
-                         PREDS1,PRED2THS3, PRED2RS3,           &
-                         PBLL_O_E,                             &
-                         PETHETA, PEMOIST                      )
-!
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
-INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
-INTEGER,                INTENT(IN)   :: KRRI          ! number of ice var.
-!
-LOGICAL,                INTENT(IN)   ::  OCLOSE_OUT   ! switch for syncronous
-                                                      ! file opening
-LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                 ! diagnostic fields in the syncronous FM-file
-CHARACTER*4           , INTENT(IN)   ::  HTURBDIM     ! Kind of turbulence param.
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                  ! metric coefficients
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF  ! Virtual Potential Temp.
-                                                  ! of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM ! Lv(T)/Cp/Exner at t-1 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM      ! Turbulent Mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS    ! Dissipative length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM,PTKEM! Conservative Potential 
-                                                  ! Temperature and TKE at t-1
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM      ! Mixing ratios at  t-1
-                                                  ! with PRM(:,:,:,1) = cons.
-                                                  ! mixing ratio
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM     ! Scalars at t-1      
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM
-                                  ! s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-!
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDTH1 ! Redelsperger number R_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDR1  ! Redelsperger number R_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2TH3 ! Redelsperger number R*2_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2R3  ! Redelsperger number R*2_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2THR3! Redelsperger number R*2_thq
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PREDS1   ! Redelsperger number R_sv
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2THS3! Redelsperger number R*2_thsv
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2RS3 ! Redelsperger number R*2_qsv
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PBLL_O_E! beta*Lk*Leps/tke
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PETHETA ! coefficient E_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PEMOIST ! coefficient E_moist
-!
-END SUBROUTINE PRANDTL
-!
-END INTERFACE
-!
-END MODULE MODI_PRANDTL
diff --git a/src/arome/turb/modi_rmc01.F90 b/src/arome/turb/modi_rmc01.F90
deleted file mode 100644
index 1845c21be2008c32a085aee0d8d689a094163dd2..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_rmc01.F90
+++ /dev/null
@@ -1,24 +0,0 @@
-!     ######spl
-      MODULE MODI_RMC01
-!     ################
-INTERFACE
-      SUBROUTINE RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY,PDZZ,PDIRCOSZW, &
-                       PSBL_DEPTH, PLMO, PLK, PLEPS                )
-!
-CHARACTER(LEN=4),         INTENT(IN)    :: HTURBLEN ! type of mixing length
-INTEGER,                  INTENT(IN)    :: KKA           !near ground array index  
-INTEGER,                  INTENT(IN)    :: KKU           !uppest atmosphere array index
-INTEGER,                  INTENT(IN)    :: KKL           !vert. levels type 1=MNH -1=ARO
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ   ! altitude of flux points
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PDXX  ! width of grid mesh (X dir)
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PDYY  ! width of grid mesh (Y dir)
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PDZZ  ! width of vert. layers
-REAL, DIMENSION(:,:),     INTENT(IN)    :: PDIRCOSZW ! Director Cosinus 
-REAL, DIMENSION(:,:),     INTENT(IN)    :: PSBL_DEPTH! SBL depth
-REAL, DIMENSION(:,:),     INTENT(IN)    :: PLMO  ! Monin Obuhkov length
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PLK   ! Mixing length
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PLEPS ! Dissipative length
-
-END SUBROUTINE RMC01
-END INTERFACE
-END MODULE MODI_RMC01
diff --git a/src/arome/turb/modi_sbl_depth.F90 b/src/arome/turb/modi_sbl_depth.F90
deleted file mode 100644
index 1f9fb94399496182fc754321a72f27ad3c047820..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_sbl_depth.F90
+++ /dev/null
@@ -1,24 +0,0 @@
-!     ######spl
-     MODULE MODI_SBL_DEPTH  
-!    ################ 
-!
-INTERFACE
-!
-      SUBROUTINE SBL_DEPTH(KKB,KKE,PZZ,PFLXU,PFLXV,PWTHV,PLMO,PSBL_DEPTH)
-!
-INTEGER,                INTENT(IN)    :: KKB       ! first physical level
-INTEGER,                INTENT(IN)    :: KKE       ! upper physical level
-REAL, DIMENSION(:,:,:), INTENT(IN)    :: PZZ       ! altitude of flux levels
-REAL, DIMENSION(:,:,:), INTENT(IN)    :: PFLXU     ! u'w'
-REAL, DIMENSION(:,:,:), INTENT(IN)    :: PFLXV     ! v'w'
-REAL, DIMENSION(:,:,:), INTENT(IN)    :: PWTHV     ! buoyancy flux
-REAL, DIMENSION(:,:),   INTENT(IN)    :: PLMO      ! Monin-Obukhov length
-REAL, DIMENSION(:,:),   INTENT(INOUT) :: PSBL_DEPTH! boundary layer height
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE SBL_DEPTH
-!
-END INTERFACE
-!
-END MODULE MODI_SBL_DEPTH
diff --git a/src/arome/turb/modi_th_r_from_thl_rt_1d.F90 b/src/arome/turb/modi_th_r_from_thl_rt_1d.F90
deleted file mode 100644
index 577edb249bc26074fe0e702f3b42755230e6f855..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_th_r_from_thl_rt_1d.F90
+++ /dev/null
@@ -1,24 +0,0 @@
-!     ######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*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
diff --git a/src/arome/turb/modi_th_r_from_thl_rt_2d.F90 b/src/arome/turb/modi_th_r_from_thl_rt_2d.F90
deleted file mode 100644
index a3af56f230d40c878f0107fccf848bf0e0767939..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_th_r_from_thl_rt_2d.F90
+++ /dev/null
@@ -1,24 +0,0 @@
-!     ######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*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
diff --git a/src/arome/turb/modi_th_r_from_thl_rt_3d.F90 b/src/arome/turb/modi_th_r_from_thl_rt_3d.F90
deleted file mode 100644
index 5b52a056f4b7b7b6889c6f0c893843bbb65cc2d7..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_th_r_from_thl_rt_3d.F90
+++ /dev/null
@@ -1,27 +0,0 @@
-!     ######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*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
diff --git a/src/arome/turb/modi_thl_rt_from_th_r_mf.F90 b/src/arome/turb/modi_thl_rt_from_th_r_mf.F90
deleted file mode 100644
index ed1dd9c6bdde9ebd379a8dbac1f5ed49281a884e..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_thl_rt_from_th_r_mf.F90
+++ /dev/null
@@ -1,31 +0,0 @@
-!     ######spl
-     MODULE MODI_THL_RT_FROM_TH_R_MF
-!    ###############################
-!
-INTERFACE
-!     #################################################################
-      SUBROUTINE THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                  &
-                                      PTH, PR, PEXN, &
-                                      PTHL, PRT                      )
-!     #################################################################
-!
-!               
-!*               1.1  Declaration of Arguments
-!                
-!
-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.
-
-REAL, DIMENSION(:,:), INTENT(IN)   :: PTH      ! theta
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PR       ! water species
-REAL, DIMENSION(:,:), INTENT(IN)   :: PEXN    ! exner function
-
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PTHL     ! th_l
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PRT      ! total non precip. water
-!
-END SUBROUTINE THL_RT_FROM_TH_R_MF
-
-END INTERFACE
-!
-END MODULE MODI_THL_RT_FROM_TH_R_MF
diff --git a/src/arome/turb/prandtl.F90 b/src/arome/turb/prandtl.F90
deleted file mode 100644
index 3afb18cc4446837c32394876d804095d204a07ea..0000000000000000000000000000000000000000
--- a/src/arome/turb/prandtl.F90
+++ /dev/null
@@ -1,501 +0,0 @@
-!     ######spl
-      SUBROUTINE PRANDTL(KKA,KKU,KKL,KRR,KRRI,OCLOSE_OUT,OTURB_DIAG,       &
-                         HTURBDIM,                             &
-                         HFMFILE,HLUOUT,                       &
-                         PDXX,PDYY,PDZZ,PDZX,PDZY,             &
-                         PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
-                         PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
-                         PREDTH1,PREDR1,                       &
-                         PRED2TH3, PRED2R3, PRED2THR3,         &
-                         PREDS1,PRED2THS3, PRED2RS3,           &
-                         PBLL_O_E,                             &
-                         PETHETA, PEMOIST                      )
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
-!     ###########################################################
-!
-!
-!!****  *PRANDTL* - routine to compute the Prandtl turbulent numbers
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this routine is to compute the Redelsperger 
-!     numbers and then get the turbulent Prandtl and Schmidt numbers:
-!       * for the heat fluxes     - PHI3 = 1/ Prandtl
-!       * for the moisture fluxes - PSI3 = 1/ Schmidt
-!
-!!**  METHOD
-!!    ------
-!!    The following steps are performed:
-!!
-!!     1 - default values of 1 are taken for phi3 and psi3 and different masks
-!!       are defined depending on the presence of turbulence, stratification and
-!!       humidity. The 1D Redelsperger numbers are computed  
-!!         * ZREDTH1 : (g / THVREF ) (LT**2 / TKE ) ETHETA (D Theta / Dz)   
-!!         * ZREDR1  : (g / THVREF ) (LT**2 / TKE ) EMOIST (D TW    / Dz)  
-!!     2 - 3D Redelsperger numbers are computed only for turbulent 
-!!       grid points where  ZREDTH1 or ZREDR1 are > 0.
-!!     3 - PHI3 is  computed only for turbulent grid points where ZREDTH1 > 0
-!!      (turbulent thermally stratified points)
-!!     4 - PSI3 is computed only for turbulent grid points where ZREDR1 > 0
-!!      (turbulent moist points)
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      FUNCTIONs ETHETA and EMOIST  :  
-!!            allows to compute the coefficients 
-!!            for the turbulent correlation between any variable
-!!            and the virtual potential temperature, of its correlations 
-!!            with the conservative potential temperature and the humidity 
-!!            conservative variable:
-!!            -------              -------              -------
-!!            A' Thv'  =  ETHETA   A' Thl'  +  EMOIST   A' Rnp'  
-!!
-!!      GX_M_M, GY_M_M, GZ_M_M :  Cartesian gradient operators
-!!      MZM : Shuman function (mean operator in the z direction)
-!!      Module MODI_ETHETA    : interface module for ETHETA
-!!      Module MODI_EMOIST    : interface module for EMOIST
-!!      Module MODI_SHUMAN    : interface module for Shuman operators
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      Module MODD_CST : contains physical constants
-!!           XG         : gravity constant
-!!
-!!      Module MODD_CTURB: contains the set of constants for
-!!                        the turbulence scheme
-!!       XCTV,XCPR2    : constants for the turbulent prandtl numbers
-!!       XTKEMIN        : minimum value allowed for the TKE
-!!
-!!      Module MODD_PARAMETERS 
-!!       JPVEXT_TURB         : number of vertical marginal points
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book 2 of documentation (routine PRANDTL)
-!!      Book 1 of documentation (Chapter: Turbulence)
-!!
-!!    AUTHOR
-!!    ------
-!!      Joan Cuxart             * INM and Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         18/10/94
-!!      Modifications: Feb 14, 1995 (J.Cuxart and J.Stein) 
-!!                                  Doctorization and Optimization
-!!      Modifications: March 21, 1995 (J.M. Carriere) 
-!!                                  Introduction of cloud water
-!!      Modifications: March 21, 1995 (J. Cuxart and J.Stein)
-!!                                  Phi3 and Psi3 at w point + cleaning
-!!      Modifications: July 2, 1995 (J.Cuxart and Ph.Bougeault)
-!!                         change the value of Phi3 and Psi3 if negative
-!!      Modifications: Sept 20, 1995 (J. Stein, J. Cuxart, J.L. Redelsperger)
-!!                         remove the Where + use REDTH1+REDR1 for the tests
-!!      Modifications: October 10, 1995 (J. Cuxart and J.Stein)
-!!                         Psi3 for tPREDS1he scalar variables
-!!      Modifications: February 27, 1996 (J.Stein) optimization
-!!      Modifications: June 15, 1996 (P.Jabouille) return to the previous
-!!                          computation of Phi3 and Psi3
-!!      Modifications: October 10, 1996 (J. Stein) change the temporal
-!!                          discretization
-!!      Modifications: May 23, 1997 (J. Stein) bug in 3D Redels number at ground
-!!                                             with orography
-!!      Modifications: Feb 20, 1998 (J. Stein) bug in all the 3D cases due to
-!!                                             the use of ZW1 instead of ZW2
-!!                     Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
-!!                     July    2005 (Tomas, Masson) implicitation of PHI3 and PSI3
-!!                     October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
-!!                                              change of YCOMMENT
-!!                     2012-02 Y. Seity,  add possibility to run with reversed 
-!!                                               vertical levels
-!!      Modifications: July 2015    (Wim de Rooy) LHARAT (Racmo turbulence) switch
-!! --------------------------------------------------------------------------
-!       
-!*      0. DECLARATIONS
-!          ------------
-!
-USE MODD_CST
-USE MODD_CONF
-USE MODD_CTURB
-USE MODD_PARAMETERS
-!
-USE MODI_GRADIENT_M
-USE MODI_EMOIST
-USE MODI_ETHETA
-USE MODI_SHUMAN, ONLY: MZM
-USE MODE_FMWRIT
-!
-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)   :: KRR           ! number of moist var.
-INTEGER,                INTENT(IN)   :: KRRI          ! number of ice var.
-!
-LOGICAL,                INTENT(IN)   ::  OCLOSE_OUT   ! switch for syncronous
-                                                      ! file opening
-LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                 ! diagnostic fields in the syncronous FM-file
-CHARACTER*4           , INTENT(IN)   ::  HTURBDIM     ! Kind of turbulence param.
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                  ! metric coefficients
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF  ! Virtual Potential Temp.
-                                                  ! of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM ! Lv(T)/Cp/Exner at t-1 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM      ! Turbulent Mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS    ! Dissipative length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM,PTKEM! Conservative Potential 
-                                                  ! Temperature and TKE at t-1
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM      ! Mixing ratios at  t-1
-                                                  ! with PRM(:,:,:,1) = cons.
-                                                  ! mixing ratio
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM     ! Scalars at t-1      
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM
-                                  ! s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-!
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDTH1 ! Redelsperger number R_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDR1  ! Redelsperger number R_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2TH3 ! Redelsperger number R*2_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2R3  ! Redelsperger number R*2_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2THR3! Redelsperger number R*2_thq
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PREDS1   ! Redelsperger number R_s
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2THS3! Redelsperger number R*2_thsv
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2RS3 ! Redelsperger number R*2_qsv
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PBLL_O_E! beta*Lk*Leps/tke
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PETHETA ! coefficient E_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PEMOIST ! coefficient E_moist
-!
-!
-!       0.2  declaration of local variables
-!
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  &
-                  ZW1, ZW2, ZW3
-!                                                 working variables
-!                                                     
-INTEGER :: IKB      ! vertical index value for the first inner mass point
-INTEGER :: IKE      ! vertical index value for the last inner mass point
-INTEGER             :: IRESP        ! Return code of FM routines
-INTEGER             :: ILENG        ! Length of the data field in LFIFM file
-INTEGER             :: IGRID        ! C-grid indicator in LFIFM file
-INTEGER             :: ILENCH       ! Length of comment string in LFIFM file
-CHARACTER (LEN=100) :: YCOMMENT     ! comment string in LFIFM file
-CHARACTER (LEN=16)  :: YRECFM       ! Name of the desired field in LFIFM file
-INTEGER::  ISV                      ! number of scalar variables       
-INTEGER::  JSV                      ! loop index for the scalar variables  
-
-INTEGER :: JLOOP
-REAL    :: ZMINVAL
-! ---------------------------------------------------------------------------
-!
-!
-!*      1.  DEFAULT VALUES,  1D REDELSPERGER NUMBERS 
-!           ----------------------------------------
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('PRANDTL',0,ZHOOK_HANDLE)
-
-IF (LHARAT) THEN
-PREDTH1(:,:,:)=0.
-PREDR1(:,:,:)=0.
-PRED2TH3(:,:,:)=0.
-PRED2R3(:,:,:)=0.
-PRED2THR3(:,:,:)=0.
-PREDS1(:,:,:,:)=0.
-PRED2THS3(:,:,:,:)=0.
-PRED2RS3(:,:,:,:)=0.
-PBLL_O_E(:,:,:)=0.
-ENDIF
-!
-IKB = KKA+JPVEXT_TURB*KKL
-IKE = KKU-JPVEXT_TURB*KKL 
-ILENG=SIZE(PTHLM,1)*SIZE(PTHLM,2)*SIZE(PTHLM,3)
-ISV  =SIZE(PSVM,4)
-!
-PETHETA(:,:,:) = MZM(ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM), KKA, KKU, KKL)
-PEMOIST(:,:,:) = MZM(EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM), KKA, KKU, KKL)
-PETHETA(:,:,KKA) = 2.*PETHETA(:,:,IKB) - PETHETA(:,:,IKB+KKL)
-PEMOIST(:,:,KKA) = 2.*PEMOIST(:,:,IKB) - PEMOIST(:,:,IKB+KKL)
-!
-!---------------------------------------------------------------------------
-IF (.NOT. LHARAT) THEN
-!
-!          1.3 1D Redelsperger numbers
-!
-PBLL_O_E(:,:,:) = MZM(XG / PTHVREF(:,:,:) * PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:), KKA, KKU, KKL)
-IF (KRR /= 0) THEN                ! moist case
-  PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * PETHETA(:,:,:) * &
-                   & GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL)
-  PREDR1(:,:,:) = XCTV*PBLL_O_E(:,:,:) * PEMOIST(:,:,:) * &
-                   & GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL)
-ELSE                              ! dry case
-  PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:)  * GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL)
-  PREDR1(:,:,:) = 0.
-END IF
-!
-!       3. Limits on 1D Redelperger numbers
-!          --------------------------------
-!
-ZMINVAL = (1.-1./XPHI_LIM)
-!
-ZW1 = 1.
-ZW2 = 1.
-!
-WHERE (PREDTH1+PREDR1<-ZMINVAL)
-  ZW1 = (-ZMINVAL) / (PREDTH1+PREDR1)
-END WHERE
-!
-WHERE (PREDTH1<-ZMINVAL)
-  ZW2 = (-ZMINVAL) / (PREDTH1)
-END WHERE
-ZW2 = MIN(ZW1,ZW2)
-!
-ZW1 = 1.
-WHERE (PREDR1<-ZMINVAL)
-  ZW1 = (-ZMINVAL) / (PREDR1)
-END WHERE
-ZW1 = MIN(ZW2,ZW1)
-!
-!
-!       3. Modification of Mixing length and dissipative length
-!          ----------------------------------------------------
-!
-PBLL_O_E(:,:,:) = PBLL_O_E(:,:,:) * ZW1(:,:,:)
-PREDTH1 (:,:,:) = PREDTH1 (:,:,:) * ZW1(:,:,:)
-PREDR1  (:,:,:) = PREDR1  (:,:,:) * ZW1(:,:,:)
-!
-!       4. Threshold for very small (in absolute value) Redelperger numbers
-!          ----------------------------------------------------------------
-!
-ZW2=SIGN(1.,PREDTH1(:,:,:))
-PREDTH1(:,:,:)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDTH1(:,:,:))
-!
-IF (KRR /= 0) THEN                ! dry case
-  ZW2=SIGN(1.,PREDR1(:,:,:))
-  PREDR1(:,:,:)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDR1(:,:,:))
-END IF
-!
-!
-!---------------------------------------------------------------------------
-!
-!          For the scalar variables
-DO JSV=1,ISV
-  PREDS1(:,:,:,JSV)=XCTV*PBLL_O_E(:,:,:)*GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
-END DO
-!
-DO JSV=1,ISV
-  ZW2=SIGN(1.,PREDS1(:,:,:,JSV))
-  PREDS1(:,:,:,JSV)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDS1(:,:,:,JSV))
-END DO
-!
-!---------------------------------------------------------------------------
-!
-!*      2.  3D REDELSPERGER NUMBERS
-!           ------------------------
-!
-IF(HTURBDIM=='1DIM') THEN        ! 1D case
-!
-!
-  PRED2TH3(:,:,:)  = PREDTH1(:,:,:)**2
-!
-  PRED2R3(:,:,:)   = PREDR1(:,:,:) **2
-!
-  PRED2THR3(:,:,:) = PREDTH1(:,:,:) * PREDR1(:,:,:)
-!
-ELSE IF (L2D) THEN                      ! 3D case in a 2D model
-!
-  IF (KRR /= 0) THEN                 ! moist 3D case
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2+(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * &
-        MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL)
-!
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-                  PEMOIST(:,:,:) * PETHETA(:,:,:) *                         &
-      MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*     &
-                     GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL), KKA, KKU, KKL)
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+KKL)
-!
-  ELSE                 ! dry 3D case in a 2D model
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  XCTV**2*PBLL_O_E(:,:,:)**2 *     &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:) = 0.
-!
-    PRED2THR3(:,:,:) = 0.
-!
-  END IF
-!
-ELSE                                 ! 3D case in a 3D model
-!
-  IF (KRR /= 0) THEN                 ! moist 3D case
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2 +  ( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 &
-      + GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 *      &
-        MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 + &
-        GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL)
-!
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-         PEMOIST(:,:,:) * PETHETA(:,:,:) *                            &
-         MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*   &
-         GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)+                           &
-         GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*                    &
-         GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL), KKA, KKU, KKL)
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+KKL)
-!
-  ELSE                 ! dry 3D case in a 3D model
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  XCTV**2*PBLL_O_E(:,:,:)**2 *                &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 &
-      + GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:) = 0.
-!
-    PRED2THR3(:,:,:) = 0.
-!
-  END IF
-!
-END IF   ! end of the if structure on the turbulence dimensionnality
-!
-!
-!---------------------------------------------------------------------------
-!
-!           5. Prandtl numbers for scalars
-!              ---------------------------
-DO JSV=1,ISV
-!
-  IF(HTURBDIM=='1DIM') THEN
-!        1D case
-    PRED2THS3(:,:,:,JSV)  = PREDS1(:,:,:,JSV) * PREDTH1(:,:,:)
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV)   = PREDR1(:,:,:) *PREDS1(:,:,:,JSV)
-    ELSE
-      PRED2RS3(:,:,:,JSV)   = 0.
-    END IF
-!
-  ELSE  IF (L2D) THEN ! 3D case in a 2D model
-!
-    IF (KRR /= 0) THEN
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL) *PETHETA
-    ELSE
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL)
-    END IF
-    PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1*                                              &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL),                 &
-                           KKA, KKU, KKL)
-!
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1 * PEMOIST *                                   &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL),          &
-                           KKA, KKU, KKL)
-    ELSE
-      PRED2RS3(:,:,:,JSV) = 0.
-    END IF
-!
-  ELSE ! 3D case in a 3D model
-!
-    IF (KRR /= 0) THEN
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL) *PETHETA
-    ELSE
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL)
-    END IF
-    PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1*                                              &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)                  &
-                          +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*       &
-                           GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL),                 &
-                           KKA, KKU, KKL)
-!
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1 * PEMOIST *                                   &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)           &
-                          +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*       &
-                           GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL),          &
-                           KKA, KKU, KKL)
-    ELSE
-      PRED2RS3(:,:,:,JSV) = 0.
-    END IF
-!
-  END IF ! end of HTURBDIM if-block
-!
-END DO
-!
-!---------------------------------------------------------------------------
-!
-!*          6. SAVES THE REDELSPERGER NUMBERS
-!              ------------------------------
-!
-IF ( OTURB_DIAG .AND. OCLOSE_OUT ) THEN
-  !
-  ! stores the RED_TH1
-  YRECFM  ='RED_TH1'
-  YCOMMENT='X_Y_Z_RED_TH1 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',PREDTH1,IGRID,ILENCH,YCOMMENT,IRESP)
-  !
-  ! stores the RED_R1
-  YRECFM  ='RED_R1'
-  YCOMMENT='X_Y_Z_RED_R1 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',PREDR1,IGRID,ILENCH,YCOMMENT,IRESP)
-  !
-  ! stores the RED2_TH3
-  YRECFM  ='RED2_TH3'
-  YCOMMENT='X_Y_Z_RED2_TH3 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',PRED2TH3,IGRID,ILENCH,YCOMMENT,IRESP)
-  !
-  ! stores the RED2_R3
-  YRECFM  ='RED2_R3'
-  YCOMMENT='X_Y_Z_RED2_R3 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',PRED2R3,IGRID,ILENCH,YCOMMENT,IRESP)
-  !
-  ! stores the RED2_THR3
-  YRECFM  ='RED2_THR3'
-  YCOMMENT='X_Y_Z_RED2_THR3 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',PRED2THR3,IGRID,ILENCH,YCOMMENT,IRESP)
-END IF
-!
-!---------------------------------------------------------------------------
-ENDIF ! (Done only if LHARAT is FALSE)
-!
-IF (LHOOK) CALL DR_HOOK('PRANDTL',1,ZHOOK_HANDLE)
-END SUBROUTINE PRANDTL
diff --git a/src/arome/turb/shallow_mf.F90 b/src/arome/turb/shallow_mf.F90
index d73839a6b6f077934c42780c6542bd5bb1d0d196..6aebb5ea30e0613b60e23d5f560dc15dbbe47013 100644
--- a/src/arome/turb/shallow_mf.F90
+++ b/src/arome/turb/shallow_mf.F90
@@ -63,7 +63,7 @@ USE MODD_CST
 USE MODD_PARAMETERS, ONLY: JPVEXT
 USE MODD_CMFSHALL
 
-USE MODI_THL_RT_FROM_TH_R_MF
+USE MODE_THL_RT_FROM_TH_R_MF, ONLY: THL_RT_FROM_TH_R_MF
 USE MODI_COMPUTE_UPDRAFT
 USE MODI_COMPUTE_UPDRAFT_RHCJ10
 USE MODI_COMPUTE_UPDRAFT_RAHA
diff --git a/src/common/turb/bl89.F90 b/src/common/turb/mode_bl89.F90
similarity index 97%
rename from src/common/turb/bl89.F90
rename to src/common/turb/mode_bl89.F90
index 57056e4815f6fb65ae1ff202472d58ed1f7810c6..9208fcb6b61d59f60206b2fac7cbf4d40ba926e0 100644
--- a/src/common/turb/bl89.F90
+++ b/src/common/turb/mode_bl89.F90
@@ -1,3 +1,10 @@
+!MNH_LIC Copyright 1994-2022 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 MODE_BL89
+IMPLICIT NONE
+CONTAINS
 !     ######spl
       SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
       USE PARKIND1, ONLY : JPRB
@@ -348,3 +355,4 @@ END IF
 !
 IF (LHOOK) CALL DR_HOOK('BL89',1,ZHOOK_HANDLE)
 END SUBROUTINE BL89
+END MODULE MODE_BL89
diff --git a/src/arome/turb/bl_depth_diag_3d.F90 b/src/common/turb/mode_bl_depth_diag.F90
similarity index 57%
rename from src/arome/turb/bl_depth_diag_3d.F90
rename to src/common/turb/mode_bl_depth_diag.F90
index 78ce7c72a1dcca18bab4db48b5a002c1b7103dfa..3cf56530a1969c0ebc35cf56bb74b7a1bc59eb1d 100644
--- a/src/arome/turb/bl_depth_diag_3d.F90
+++ b/src/common/turb/mode_bl_depth_diag.F90
@@ -1,4 +1,16 @@
-!     ######spl
+!MNH_LIC Copyright 1994-2022 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 MODE_BL_DEPTH_DIAG
+!
+INTERFACE BL_DEPTH_DIAG  
+      MODULE PROCEDURE BL_DEPTH_DIAG_3D
+      MODULE PROCEDURE BL_DEPTH_DIAG_1D
+END INTERFACE
+!
+CONTAINS
+!
 FUNCTION BL_DEPTH_DIAG_3D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -91,3 +103,41 @@ BL_DEPTH_DIAG_3D(:,:) = BL_DEPTH_DIAG_3D(:,:) / (1. - PFTOP_O_FSURF)
 !
 IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_3D',1,ZHOOK_HANDLE)
 END FUNCTION BL_DEPTH_DIAG_3D
+!
+FUNCTION BL_DEPTH_DIAG_1D(KKB,KKE,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF)
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+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
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',0,ZHOOK_HANDLE)
+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)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',1,ZHOOK_HANDLE)
+END FUNCTION BL_DEPTH_DIAG_1D
+END MODULE MODE_BL_DEPTH_DIAG
diff --git a/src/common/turb/emoist.F90 b/src/common/turb/mode_emoist.F90
similarity index 97%
rename from src/common/turb/emoist.F90
rename to src/common/turb/mode_emoist.F90
index f17bc3168c165b43e3bac222dad961b974e4834f..b06579b1f695af14ad4dc9232fcaab9e9b456bee 100644
--- a/src/common/turb/emoist.F90
+++ b/src/common/turb/mode_emoist.F90
@@ -1,8 +1,10 @@
-!MNH_LIC Copyright 1995-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-2022 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 MODE_EMOIST
+IMPLICIT NONE
+CONTAINS
 FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -162,3 +164,4 @@ END IF
 !
 IF (LHOOK) CALL DR_HOOK('EMOIST',1,ZHOOK_HANDLE)
 END FUNCTION EMOIST
+END MODULE MODE_EMOIST
diff --git a/src/common/turb/etheta.F90 b/src/common/turb/mode_etheta.F90
similarity index 98%
rename from src/common/turb/etheta.F90
rename to src/common/turb/mode_etheta.F90
index f0506bd89dc5fc455e55281835ba7c1ab6b2365c..4e3e91fad378230da49f3c590fc0784a50b16144 100644
--- a/src/common/turb/etheta.F90
+++ b/src/common/turb/mode_etheta.F90
@@ -2,7 +2,9 @@
 !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 MODE_ETHETA
+IMPLICIT NONE
+CONTAINS
 FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -153,3 +155,4 @@ END IF
 !
 IF (LHOOK) CALL DR_HOOK('ETHETA',1,ZHOOK_HANDLE)
 END FUNCTION ETHETA
+END MODULE MODE_ETHETA
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index 6a84eb9b3cb90b711360c9ffa7f765e65ca905b8..9cce0c8742d08deed5e47cc3254740f5eafbde8e 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -151,8 +151,8 @@ USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 !
 USE MODI_GRADIENT_M
-USE MODI_EMOIST
-USE MODI_ETHETA
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
 USE MODI_SHUMAN, ONLY: MZM
 USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 !
diff --git a/src/common/turb/rmc01.F90 b/src/common/turb/mode_rmc01.F90
similarity index 95%
rename from src/common/turb/rmc01.F90
rename to src/common/turb/mode_rmc01.F90
index 7488d6cbd8f5651433d9f34838940920a54ea5f4..628b4cad0dda1fcef65f2a70fd43b0ef5ec61e0f 100644
--- a/src/common/turb/rmc01.F90
+++ b/src/common/turb/mode_rmc01.F90
@@ -1,5 +1,11 @@
-!     ######spl
-      SUBROUTINE RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY, &
+!MNH_LIC Copyright 1994-2022 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 MODE_RMC01
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY, &
       PDZZ,PDIRCOSZW,PSBL_DEPTH,PLMO,PLK,PLEPS)
       USE PARKIND1, ONLY : JPRB
       USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -233,3 +239,4 @@ PLEPS(:,:,KKU  ) = PLEPS(:,:,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('RMC01',1,ZHOOK_HANDLE)
 END SUBROUTINE RMC01
+END MODULE MODE_RMC01
diff --git a/src/common/turb/sbl_depth.F90 b/src/common/turb/mode_sbl_depth.F90
similarity index 91%
rename from src/common/turb/sbl_depth.F90
rename to src/common/turb/mode_sbl_depth.F90
index 0c670f8db17dc556afeb5c44ab5df47793acf41d..e485940e06f075776b2f2ddd41aac71164aa331b 100644
--- a/src/common/turb/sbl_depth.F90
+++ b/src/common/turb/mode_sbl_depth.F90
@@ -1,3 +1,10 @@
+!MNH_LIC Copyright 1994-2022 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 MODE_SBL_DEPTH
+IMPLICIT NONE
+CONTAINS
 !     ######spl
       SUBROUTINE SBL_DEPTH(KKB,KKE,PZZ,PFLXU,PFLXV,PWTHV,PLMO,PSBL_DEPTH)
       USE PARKIND1, ONLY : JPRB
@@ -41,7 +48,7 @@
 USE MODD_PARAMETERS, ONLY : XUNDEF
 USE MODD_CTURB,      ONLY : XFTOP_O_FSURF, XSBL_O_BL
 !
-USE MODI_BL_DEPTH_DIAG
+USE MODE_BL_DEPTH_DIAG
 !
 IMPLICIT NONE
 !
@@ -119,3 +126,4 @@ WHERE (PLMO(:,:)==XUNDEF) PSBL_DEPTH = ZSBL_DYN
 !----------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('SBL_DEPTH',1,ZHOOK_HANDLE)
 END SUBROUTINE SBL_DEPTH
+END MODULE MODE_SBL_DEPTH
diff --git a/src/arome/turb/th_r_from_thl_rt_1d.F90 b/src/common/turb/mode_th_r_from_thl_rt_1d.F90
similarity index 92%
rename from src/arome/turb/th_r_from_thl_rt_1d.F90
rename to src/common/turb/mode_th_r_from_thl_rt_1d.F90
index f553f45851e617143f842eca05a70f4eb9817456..d042e376d5a6cfa46f04793cf75bc0908726113a 100644
--- a/src/arome/turb/th_r_from_thl_rt_1d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_1d.F90
@@ -1,4 +1,10 @@
-!     ######spl
+!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_1D
+IMPLICIT NONE
+CONTAINS
       SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
                                   PTHL, PRT, PTH, PRV, PRL, PRI,         &
                                   PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
@@ -35,6 +41,7 @@
 !!      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'
 !!
 !! --------------------------------------------------------------------------
 !
@@ -45,6 +52,7 @@ USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 USE MODE_COMPUTE_FRAC_ICE, ONLY : COMPUTE_FRAC_ICE
 USE MODD_CST!, ONLY: XP00, XRD, XCPD, XCPV, XCL, XCI, XLVTT, XTT, XLSTT
+USE MODD_DYN_n, ONLY : LOCEAN
 USE MODE_THERMO
 !
 IMPLICIT NONE
@@ -52,7 +60,7 @@ IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
 !
-CHARACTER*1         , INTENT(IN) :: HFRAC_ICE
+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
@@ -119,8 +127,11 @@ ENDDO
 !         ---------
 
 DO II=1,JITER
-  ZT(:)=PTH(:)*ZEXN(:)
-
+  IF (LOCEAN) THEN
+    ZT=PTH                  
+  ELSE
+    ZT(:)=PTH(:)*ZEXN(:)
+  END IF
   !Computation of liquid/ice fractions
   PFRAC_ICE(:) = 0.
   DO J=1, SIZE(PFRAC_ICE, 1)
@@ -192,3 +203,4 @@ ENDDO
 IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_1D',1,ZHOOK_HANDLE)
 
 END SUBROUTINE TH_R_FROM_THL_RT_1D
+END MODULE MODE_TH_R_FROM_THL_RT_1D
diff --git a/src/arome/turb/th_r_from_thl_rt_2d.F90 b/src/common/turb/mode_th_r_from_thl_rt_2d.F90
similarity index 87%
rename from src/arome/turb/th_r_from_thl_rt_2d.F90
rename to src/common/turb/mode_th_r_from_thl_rt_2d.F90
index 5d1ff0e0802183504784a6b10abb3c76943c72df..e4292fa121fc7d17360a06e5e625b7630f978487 100644
--- a/src/arome/turb/th_r_from_thl_rt_2d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_2d.F90
@@ -1,4 +1,10 @@
-!     ######spl
+!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_2D
+IMPLICIT NONE
+CONTAINS
       SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP,             &
                                   PTHL, PRT, PTH, PRV, PRL, PRI,         &
                                   PRSATW, PRSATI, PRR, PRS, PRG, PRH     )
@@ -40,7 +46,7 @@
 !          ------------
 !
 !
-USE MODI_TH_R_FROM_THL_RT_3D
+USE MODE_TH_R_FROM_THL_RT_3D, ONLY: TH_R_FROM_THL_RT_3D
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 
@@ -49,7 +55,7 @@ IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
 !
-CHARACTER*1           , INTENT(IN) :: HFRAC_ICE
+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.
@@ -101,3 +107,5 @@ ENDDO
 IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_2D',1,ZHOOK_HANDLE)
 
 END SUBROUTINE TH_R_FROM_THL_RT_2D
+END MODULE MODE_TH_R_FROM_THL_RT_2D
+
diff --git a/src/arome/turb/th_r_from_thl_rt_3d.F90 b/src/common/turb/mode_th_r_from_thl_rt_3d.F90
similarity index 87%
rename from src/arome/turb/th_r_from_thl_rt_3d.F90
rename to src/common/turb/mode_th_r_from_thl_rt_3d.F90
index 473a9bcc278263fd8ea78fed112a4520567b0f21..fb42b71011f9b8c3029ad443ec2f3a036ec6caad 100644
--- a/src/arome/turb/th_r_from_thl_rt_3d.F90
+++ b/src/common/turb/mode_th_r_from_thl_rt_3d.F90
@@ -1,4 +1,10 @@
-!     ######spl
+!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_3D
+IMPLICIT NONE
+CONTAINS      
       SUBROUTINE TH_R_FROM_THL_RT_3D(HFRAC_ICE,PFRAC_ICE,PP,             &
                                   PTHL, PRT, PTH, PRV, PRL, PRI, &
                                   PRSATW, PRSATI, PRR, PRS, PRG, PRH       )
@@ -39,7 +45,7 @@
 !*      0. DECLARATIONS
 !          ------------
 !
-USE MODI_TH_R_FROM_THL_RT_1D
+USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !
@@ -48,7 +54,7 @@ IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
 !
-CHARACTER*1         , INTENT(IN) :: HFRAC_ICE
+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
@@ -98,3 +104,4 @@ ENDDO
 IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_3D',1,ZHOOK_HANDLE)
 
 END SUBROUTINE TH_R_FROM_THL_RT_3D
+END MODULE MODE_TH_R_FROM_THL_RT_3D
diff --git a/src/arome/turb/thl_rt_from_th_r_mf.F90 b/src/common/turb/mode_thl_rt_from_th_r_mf.F90
similarity index 90%
rename from src/arome/turb/thl_rt_from_th_r_mf.F90
rename to src/common/turb/mode_thl_rt_from_th_r_mf.F90
index c9488f86b61a8c1e520804f2a283ccd26e28a208..bf72ab9439b62e883471b54f7abb10ecce8c4031 100644
--- a/src/arome/turb/thl_rt_from_th_r_mf.F90
+++ b/src/common/turb/mode_thl_rt_from_th_r_mf.F90
@@ -1,4 +1,10 @@
-!     ######spl
+!MNH_LIC Copyright 1994-2022 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 MODE_THL_RT_FROM_TH_R_MF
+IMPLICIT NONE
+CONTAINS
       SUBROUTINE THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                  &
                                       PTH, PR, PEXN, &
                                       PTHL, PRT                      )
@@ -114,3 +120,4 @@ ELSE
 END IF
 IF (LHOOK) CALL DR_HOOK('THL_RT_FRM_TH_R_MF',1,ZHOOK_HANDLE)
 END SUBROUTINE THL_RT_FROM_TH_R_MF
+END MODULE MODE_THL_RT_FROM_TH_R_MF
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index 2f31473d774a9b4cb9aa380c8a67e63cd7405acd..a55c661ce21929b9f0705255e9b821303f6a8a65 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -241,11 +241,6 @@ IKE=KKU-JPVEXT_TURB*KKL
 ! compute the effective diffusion coefficient at the mass point
 ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) 
 !
-#ifdef REPRO48
-#else
-IF (LBUDGET_TH)  CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH),  'DISSH', PRTHLS(:, :, :) )
-#endif
-!
 !----------------------------------------------------------------------------
 !
 !*       2.   TKE EQUATION  
@@ -383,11 +378,6 @@ IF (LBUDGET_TKE) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TKE), 'TR', PRTKES(:, :
 !
 PRTHLS(:,:,:) = PRTHLS(:,:,:) + XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * &
                 (PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:) * PCOEF_DISS(:,:,:)
-
-#ifdef REPRO48
-#else
-IF (LBUDGET_TH) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
-#endif
 !----------------------------------------------------------------------------
 !
 !*       4.   STORES SOME DIAGNOSTICS
diff --git a/src/common/turb/mode_turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
index 32e00a0afc6051fc5a73e7608a4773a969e20e4a..2860adfd30c427a8fe410d348ff76ca29733595a 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -219,8 +219,8 @@ 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
@@ -230,7 +230,7 @@ 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/common/turb/mode_turb_ver_sv_corr.F90 b/src/common/turb/mode_turb_ver_sv_corr.F90
index 7de442316cf8bdb5210ef16ce366845d4da4eb8a..97928bf69e4f3fc0928ccc66342956685e3a14c2 100644
--- a/src/common/turb/mode_turb_ver_sv_corr.F90
+++ b/src/common/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/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index 6fdbba5c3f46700856f79ee23f60923d99b19c5e..781e1b125b46cd8c16d3dd9311d7803201ac0fd8 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/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/arome/turb/modi_turb.F90 b/src/common/turb/modi_turb.F90
similarity index 79%
rename from src/arome/turb/modi_turb.F90
rename to src/common/turb/modi_turb.F90
index 0120f583bdfdb428e3fe41b84f62438a765d36af..032c257502511d07f5ae84ce9133d9b93d9c9239 100644
--- a/src/arome/turb/modi_turb.F90
+++ b/src/common/turb/modi_turb.F90
@@ -6,18 +6,18 @@ INTERFACE
 !
       SUBROUTINE TURB(KKA, KKU, KKL, KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,   &
               & KSPLIT,KMODEL_CL, &
-              & OCLOSE_OUT,OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
+              & OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
               & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HINST_SFU,         &
-              & HMF_UPDRAFT,PIMPL,PTSTEP_UVW, PTSTEP_MET,PTSTEP_SV,   &
-              & HFMFILE,HLUOUT,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,          &
+              & HMF_UPDRAFT,PIMPL,PTSTEP,TPFILE,                      &
+              & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,          &
               & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
               & PRHODJ,PTHVREF,PRHODREF,                              &
               & PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
-              & PPABSM,PUM,PVM,PWM,PTKEM,PSVM,PSRCM,                  &
+              & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
               & PLENGTHM,PLENGTHH,MFMOIST,                            &
               & PBL_DEPTH, PSBL_DEPTH,                                &
-              & PUT,PVT,PWT,PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,    &
-              & PTHLM,PRM,                                            &
+              & PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,    &
+              & PTHLT,PRT,                                            &
               & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,              &
               & PHGRAD,PSIGS,                                         &
               & PDRUS_TURB,PDRVS_TURB,                                &
@@ -25,12 +25,13 @@ INTERFACE
               & PFLXZTHVMF,PWTH,PWRC,PWSV,PDP,PTP,PTPMF,PTDIFF,PTDISS,&
               & YDDDH,YDLDDH,YDMDDH,                                  &
               & TBUDGETS, KBUDGETS,                                   &
-              & PTR,PDISS,PEDR                                        ) 
+              & PTR,PDISS,PEDR,PLEM                                   ) 
 !
 USE DDH_MIX, ONLY  : TYP_DDH
 USE YOMLDDH, ONLY  : TLDDH
 USE YOMMDDH, ONLY  : TMDDH
 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
@@ -43,8 +44,6 @@ CHARACTER(LEN=*),DIMENSION(:),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
 CHARACTER(LEN=4),INTENT(IN)          :: HMF_UPDRAFT   ! Type of mass flux
 INTEGER,                INTENT(IN)   :: KSPLIT        ! number of time-splitting
 INTEGER,                INTENT(IN)   :: KMODEL_CL     ! model number for cloud mixing length
-LOGICAL,                INTENT(IN)   ::  OCLOSE_OUT   ! switch for syncronous
-                                                      ! file opening
 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
@@ -52,21 +51,16 @@ LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
 LOGICAL,                INTENT(IN)   ::  OSUBG_COND   ! switch for SUBGrid 
                                  ! CONDensation
 LOGICAL,                INTENT(IN)   ::  ORMC01       ! switch for RMC01 lengths in SBL
-CHARACTER*4           , INTENT(IN)      ::  HTURBDIM  ! dimensionality of the 
+CHARACTER(LEN=4)       , INTENT(IN)   ::  HTURBDIM  ! dimensionality of the 
                                  ! turbulence scheme
-CHARACTER*4           , INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
-CHARACTER*4           , INTENT(IN)   ::  HTOM         ! kind of Third Order Moment
-CHARACTER*4           , INTENT(IN)   ::  HTURBLEN_CL  ! kind of cloud mixing length
-CHARACTER*1           , INTENT(IN)   ::  HINST_SFU    ! temporal location of the
+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
+CHARACTER(LEN=1)      , INTENT(IN)   ::  HINST_SFU    ! temporal location of the
                                                       ! surface friction flux
 REAL,                   INTENT(IN)   ::  PIMPL        ! degree of implicitness
-REAL,                   INTENT(IN)   ::  PTSTEP_UVW   ! Dynamical timestep 
-REAL,                   INTENT(IN)   ::  PTSTEP_MET   ! Timestep for meteorological variables                        
-REAL,                   INTENT(IN)   ::  PTSTEP_SV    ! Timestep for tracer variables
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+REAL,                   INTENT(IN)   ::  PTSTEP       ! Timestep
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file for MesoNH
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDXX,PDYY,PDZZ,PDZX,PDZY
                                         ! metric coefficients
@@ -94,17 +88,15 @@ REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PSFSV
 ! normal surface fluxes of Scalar var. 
 !
 !    prognostic variables at t- deltat
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM,PWM ! wind components
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM        ! passive scal. var.
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PSRCM       ! Second-order flux
+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
 !
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PUT,PVT,PWT ! Wind  at t
-!
 !    variables for cloud mixing length
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PCEI ! Cloud Entrainment instability
                                                  ! index to emphasize localy 
@@ -113,8 +105,8 @@ REAL, INTENT(IN)      ::  PCEI_MIN ! minimum threshold for the instability index
 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) ::  PTHLM       ! conservative pot. temp.
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRM         ! water var.  where 
+REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PTHLT       ! conservative pot. temp.
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRT         ! water var.  where 
                              ! PRM(:,:,:,1) is the conservative mixing ratio        
 !
 ! sources of momentum, conservative potential temperature, Turb. Kin. Energy, 
@@ -158,7 +150,8 @@ INTEGER, INTENT(IN) :: KBUDGETS
 !
 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  :: PEDR  ! EDR
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PLEM  ! Mixing length
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/common/turb/prandtl.F90 b/src/common/turb/prandtl.F90
deleted file mode 100644
index 458899f3908ef3d12ae9a65efd25f416b8245f46..0000000000000000000000000000000000000000
--- a/src/common/turb/prandtl.F90
+++ /dev/null
@@ -1,540 +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 MODE_PRANDTL
-IMPLICIT NONE
-CONTAINS
-      SUBROUTINE PRANDTL(KKA,KKU,KKL,KRR,KRRI,OTURB_DIAG,       &
-                         HTURBDIM,                             &
-                         TPFILE,                               &
-                         PDXX,PDYY,PDZZ,PDZX,PDZY,             &
-                         PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
-                         PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
-                         PREDTH1,PREDR1,                       &
-                         PRED2TH3, PRED2R3, PRED2THR3,         &
-                         PREDS1,PRED2THS3, PRED2RS3,           &
-                         PBLL_O_E,                             &
-                         PETHETA, PEMOIST                      )
-!     ###########################################################
-!
-!
-!!****  *PRANDTL* - routine to compute the Prandtl turbulent numbers
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this routine is to compute the Redelsperger 
-!     numbers and then get the turbulent Prandtl and Schmidt numbers:
-!       * for the heat fluxes     - PHI3 = 1/ Prandtl
-!       * for the moisture fluxes - PSI3 = 1/ Schmidt
-!
-!!**  METHOD
-!!    ------
-!!    The following steps are performed:
-!!
-!!     1 - default values of 1 are taken for phi3 and psi3 and different masks
-!!       are defined depending on the presence of turbulence, stratification and
-!!       humidity. The 1D Redelsperger numbers are computed  
-!!         * ZREDTH1 : (g / THVREF ) (LT**2 / TKE ) ETHETA (D Theta / Dz)   
-!!         * ZREDR1  : (g / THVREF ) (LT**2 / TKE ) EMOIST (D TW    / Dz)  
-!!     2 - 3D Redelsperger numbers are computed only for turbulent 
-!!       grid points where  ZREDTH1 or ZREDR1 are > 0.
-!!     3 - PHI3 is  computed only for turbulent grid points where ZREDTH1 > 0
-!!      (turbulent thermally stratified points)
-!!     4 - PSI3 is computed only for turbulent grid points where ZREDR1 > 0
-!!      (turbulent moist points)
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      FUNCTIONs ETHETA and EMOIST  :  
-!!            allows to compute the coefficients 
-!!            for the turbulent correlation between any variable
-!!            and the virtual potential temperature, of its correlations 
-!!            with the conservative potential temperature and the humidity 
-!!            conservative variable:
-!!            -------              -------              -------
-!!            A' Thv'  =  ETHETA   A' Thl'  +  EMOIST   A' Rnp'  
-!!
-!!      GX_M_M, GY_M_M, GZ_M_M :  Cartesian gradient operators
-!!      MZM : Shuman function (mean operator in the z direction)
-!!      Module MODI_ETHETA    : interface module for ETHETA
-!!      Module MODI_EMOIST    : interface module for EMOIST
-!!      Module MODI_SHUMAN    : interface module for Shuman operators
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      Module MODD_CST : contains physical constants
-!!           XG         : gravity constant
-!!
-!!      Module MODD_CTURB: contains the set of constants for
-!!                        the turbulence scheme
-!!       XCTV,XCPR2    : constants for the turbulent prandtl numbers
-!!       XTKEMIN        : minimum value allowed for the TKE
-!!
-!!      Module MODD_PARAMETERS 
-!!       JPVEXT_TURB         : number of vertical marginal points
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book 2 of documentation (routine PRANDTL)
-!!      Book 1 of documentation (Chapter: Turbulence)
-!!
-!!    AUTHOR
-!!    ------
-!!      Joan Cuxart             * INM and Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original         18/10/94
-!!      Modifications: Feb 14, 1995 (J.Cuxart and J.Stein) 
-!!                                  Doctorization and Optimization
-!!      Modifications: March 21, 1995 (J.M. Carriere) 
-!!                                  Introduction of cloud water
-!!      Modifications: March 21, 1995 (J. Cuxart and J.Stein)
-!!                                  Phi3 and Psi3 at w point + cleaning
-!!      Modifications: July 2, 1995 (J.Cuxart and Ph.Bougeault)
-!!                         change the value of Phi3 and Psi3 if negative
-!!      Modifications: Sept 20, 1995 (J. Stein, J. Cuxart, J.L. Redelsperger)
-!!                         remove the Where + use REDTH1+REDR1 for the tests
-!!      Modifications: October 10, 1995 (J. Cuxart and J.Stein)
-!!                         Psi3 for tPREDS1he scalar variables
-!!      Modifications: February 27, 1996 (J.Stein) optimization
-!!      Modifications: June 15, 1996 (P.Jabouille) return to the previous
-!!                          computation of Phi3 and Psi3
-!!      Modifications: October 10, 1996 (J. Stein) change the temporal
-!!                          discretization
-!!      Modifications: May 23, 1997 (J. Stein) bug in 3D Redels number at ground
-!!                                             with orography
-!!      Modifications: Feb 20, 1998 (J. Stein) bug in all the 3D cases due to
-!!                                             the use of ZW1 instead of ZW2
-!!                     Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
-!!                     July    2005 (Tomas, Masson) implicitation of PHI3 and PSI3
-!!                     October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
-!!                                              change of YCOMMENT
-!!                     2012-02 Y. Seity,  add possibility to run with reversed 
-!!                                               vertical levels
-!!      Modifications: July 2015    (Wim de Rooy) LHARAT (Racmo turbulence) switch
-!!                     2017-09 J.Escobar, use epsilon XMNH_TINY_12 for R*4 
-!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!! JL Redelsperger 03/2021 : adding Ocean case for temperature only 
-!! --------------------------------------------------------------------------
-!       
-!*      0. DECLARATIONS
-!          ------------
-!
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-!
-USE MODD_CST
-USE MODD_CONF
-USE MODD_CTURB
-USE MODD_DYN_n,          ONLY: LOCEAN
-USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
-USE MODD_IO,             ONLY: TFILEDATA
-USE MODD_PARAMETERS
-!
-USE MODI_GRADIENT_M
-USE MODI_EMOIST
-USE MODI_ETHETA
-USE MODI_SHUMAN, ONLY: MZM
-USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
-!
-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)   :: KRR           ! number of moist var.
-INTEGER,                INTENT(IN)   :: KRRI          ! number of ice var.
-!
-LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                 ! diagnostic fields in the syncronous FM-file
-CHARACTER(LEN=4),       INTENT(IN)   ::  HTURBDIM     ! Kind of turbulence param.
-TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                  ! metric coefficients
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF  ! Virtual Potential Temp.
-                                                  ! of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM ! Lv(T)/Cp/Exner at t-1 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM      ! Turbulent Mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS    ! Dissipative length
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM,PTKEM! Conservative Potential 
-                                                  ! Temperature and TKE at t-1
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM      ! Mixing ratios at  t-1
-                                                  ! with PRM(:,:,:,1) = cons.
-                                                  ! mixing ratio
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM     ! Scalars at t-1      
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM
-                                  ! s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-!
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDTH1 ! Redelsperger number R_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PREDR1  ! Redelsperger number R_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2TH3 ! Redelsperger number R*2_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2R3  ! Redelsperger number R*2_q
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PRED2THR3! Redelsperger number R*2_thq
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PREDS1   ! Redelsperger number R_s
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2THS3! Redelsperger number R*2_thsv
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)::  PRED2RS3 ! Redelsperger number R*2_qsv
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PBLL_O_E! beta*Lk*Leps/tke
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PETHETA ! coefficient E_theta
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PEMOIST ! coefficient E_moist
-!
-!
-!       0.2  declaration of local variables
-!
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  &
-                  ZW1, ZW2, ZW3
-!                                                 working variables
-!                                                     
-INTEGER :: IKB      ! vertical index value for the first inner mass point
-INTEGER :: IKE      ! vertical index value for the last inner mass point
-INTEGER             :: IRESP        ! Return code of FM routines
-INTEGER             :: ILENG        ! Length of the data field in LFIFM file
-INTEGER             :: IGRID        ! C-grid indicator in LFIFM file
-INTEGER             :: ILENCH       ! Length of comment string in LFIFM file
-CHARACTER (LEN=100) :: YCOMMENT     ! comment string in LFIFM file
-CHARACTER (LEN=16)  :: YRECFM       ! Name of the desired field in LFIFM file
-INTEGER::  ISV                      ! number of scalar variables       
-INTEGER::  JSV                      ! loop index for the scalar variables  
-
-INTEGER :: JLOOP
-REAL    :: ZMINVAL
-TYPE(TFIELDDATA)  :: TZFIELD
-! ---------------------------------------------------------------------------
-!
-!*      1.  DEFAULT VALUES,  1D REDELSPERGER NUMBERS 
-!           ----------------------------------------
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('PRANDTL',0,ZHOOK_HANDLE)
-
-IF (LHARAT) THEN
-PREDTH1(:,:,:)=0.
-PREDR1(:,:,:)=0.
-PRED2TH3(:,:,:)=0.
-PRED2R3(:,:,:)=0.
-PRED2THR3(:,:,:)=0.
-PREDS1(:,:,:,:)=0.
-PRED2THS3(:,:,:,:)=0.
-PRED2RS3(:,:,:,:)=0.
-PBLL_O_E(:,:,:)=0.
-ENDIF
-!
-IKB = KKA+JPVEXT_TURB*KKL
-IKE = KKU-JPVEXT_TURB*KKL 
-ILENG=SIZE(PTHLM,1)*SIZE(PTHLM,2)*SIZE(PTHLM,3)
-ISV  =SIZE(PSVM,4)
-!
-PETHETA(:,:,:) = MZM(ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM), KKA, KKU, KKL)
-PEMOIST(:,:,:) = MZM(EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM), KKA, KKU, KKL)
-PETHETA(:,:,KKA) = 2.*PETHETA(:,:,IKB) - PETHETA(:,:,IKB+KKL)
-PEMOIST(:,:,KKA) = 2.*PEMOIST(:,:,IKB) - PEMOIST(:,:,IKB+KKL)
-!
-!---------------------------------------------------------------------------
-IF (.NOT. LHARAT) THEN
-!
-!          1.3 1D Redelsperger numbers
-!
-PBLL_O_E(:,:,:) = MZM(XG / PTHVREF(:,:,:) * PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:), KKA, KKU, KKL)
-IF (KRR /= 0) THEN                ! moist case
-  PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * PETHETA(:,:,:) * &
-                   & GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL)
-  PREDR1(:,:,:) = XCTV*PBLL_O_E(:,:,:) * PEMOIST(:,:,:) * &
-                   & GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL)
-ELSE                              ! dry case
-  PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:)  * GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL)
-  PREDR1(:,:,:) = 0.
-END IF
-!
-!       3. Limits on 1D Redelperger numbers
-!          --------------------------------
-!
-ZMINVAL = (1.-1./XPHI_LIM)
-!
-ZW1 = 1.
-ZW2 = 1.
-!
-WHERE (PREDTH1+PREDR1<-ZMINVAL)
-  ZW1 = (-ZMINVAL) / (PREDTH1+PREDR1)
-END WHERE
-!
-WHERE (PREDTH1<-ZMINVAL)
-  ZW2 = (-ZMINVAL) / (PREDTH1)
-END WHERE
-ZW2 = MIN(ZW1,ZW2)
-!
-ZW1 = 1.
-WHERE (PREDR1<-ZMINVAL)
-  ZW1 = (-ZMINVAL) / (PREDR1)
-END WHERE
-ZW1 = MIN(ZW2,ZW1)
-!
-!
-!       3. Modification of Mixing length and dissipative length
-!          ----------------------------------------------------
-!
-PBLL_O_E(:,:,:) = PBLL_O_E(:,:,:) * ZW1(:,:,:)
-PREDTH1 (:,:,:) = PREDTH1 (:,:,:) * ZW1(:,:,:)
-PREDR1  (:,:,:) = PREDR1  (:,:,:) * ZW1(:,:,:)
-!
-!       4. Threshold for very small (in absolute value) Redelperger numbers
-!          ----------------------------------------------------------------
-!
-ZW2=SIGN(1.,PREDTH1(:,:,:))
-PREDTH1(:,:,:)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDTH1(:,:,:))
-!
-IF (KRR /= 0) THEN                ! dry case
-  ZW2=SIGN(1.,PREDR1(:,:,:))
-  PREDR1(:,:,:)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDR1(:,:,:))
-END IF
-!
-!
-!---------------------------------------------------------------------------
-!
-!          For the scalar variables
-DO JSV=1,ISV
-  PREDS1(:,:,:,JSV)=XCTV*PBLL_O_E(:,:,:)*GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
-END DO
-!
-DO JSV=1,ISV
-  ZW2=SIGN(1.,PREDS1(:,:,:,JSV))
-  PREDS1(:,:,:,JSV)= ZW2(:,:,:) * MAX(1.E-30, ZW2(:,:,:)*PREDS1(:,:,:,JSV))
-END DO
-!
-!---------------------------------------------------------------------------
-!
-!*      2.  3D REDELSPERGER NUMBERS
-!           ------------------------
-!
-IF(HTURBDIM=='1DIM') THEN        ! 1D case
-!
-!
-  PRED2TH3(:,:,:)  = PREDTH1(:,:,:)**2
-!
-  PRED2R3(:,:,:)   = PREDR1(:,:,:) **2
-!
-  PRED2THR3(:,:,:) = PREDTH1(:,:,:) * PREDR1(:,:,:)
-!
-ELSE IF (L2D) THEN                      ! 3D case in a 2D model
-!
-  IF (KRR /= 0) THEN                 ! moist 3D case
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2+(XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * &
-        MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL)
-!
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-                  PEMOIST(:,:,:) * PETHETA(:,:,:) *                         &
-      MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*     &
-                     GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL), KKA, KKU, KKL)
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+KKL)
-!
-  ELSE                 ! dry 3D case in a 2D model
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  XCTV**2*PBLL_O_E(:,:,:)**2 *     &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:) = 0.
-!
-    PRED2THR3(:,:,:) = 0.
-!
-  END IF
-!
-ELSE                                 ! 3D case in a 3D model
-!
-  IF (KRR /= 0) THEN                 ! moist 3D case
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2 +  ( XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 &
-      + GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 *      &
-        MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 + &
-        GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+KKL)
-!
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-         PEMOIST(:,:,:) * PETHETA(:,:,:) *                            &
-         MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*   &
-         GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)+                           &
-         GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*                    &
-         GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL), KKA, KKU, KKL)
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+KKL)
-!
-  ELSE                 ! dry 3D case in a 3D model
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  XCTV**2*PBLL_O_E(:,:,:)**2 *                &
-      MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)**2 &
-      + GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)**2, KKA, KKU, KKL)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+KKL)
-!
-    PRED2R3(:,:,:) = 0.
-!
-    PRED2THR3(:,:,:) = 0.
-!
-  END IF
-!
-END IF   ! end of the if structure on the turbulence dimensionnality
-!
-!
-!---------------------------------------------------------------------------
-!
-!           5. Prandtl numbers for scalars
-!              ---------------------------
-DO JSV=1,ISV
-!
-  IF(HTURBDIM=='1DIM') THEN
-!        1D case
-    PRED2THS3(:,:,:,JSV)  = PREDS1(:,:,:,JSV) * PREDTH1(:,:,:)
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV)   = PREDR1(:,:,:) *PREDS1(:,:,:,JSV)
-    ELSE
-      PRED2RS3(:,:,:,JSV)   = 0.
-    END IF
-!
-  ELSE  IF (L2D) THEN ! 3D case in a 2D model
-!
-    IF (KRR /= 0) THEN
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL) *PETHETA
-    ELSE
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL)
-    END IF
-    PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1*                                              &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL),                 &
-                           KKA, KKU, KKL)
-!
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1 * PEMOIST *                                   &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL),          &
-                           KKA, KKU, KKL)
-    ELSE
-      PRED2RS3(:,:,:,JSV) = 0.
-    END IF
-!
-  ELSE ! 3D case in a 3D model
-!
-    IF (KRR /= 0) THEN
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL) *PETHETA
-    ELSE
-      ZW1 = MZM((XG / PTHVREF * PLM * PLEPS / PTKEM)**2, KKA, KKU, KKL)
-    END IF
-    PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1*                                              &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)                  &
-                          +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*       &
-                           GY_M_M(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL),                 &
-                           KKA, KKU, KKL)
-!
-    IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                       ZW1 * PEMOIST *                                   &
-                       MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, KKA, KKU, KKL)*       &
-                           GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)           &
-                          +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, KKA, KKU, KKL)*       &
-                           GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, KKA, KKU, KKL),          &
-                           KKA, KKU, KKL)
-    ELSE
-      PRED2RS3(:,:,:,JSV) = 0.
-    END IF
-!
-  END IF ! end of HTURBDIM if-block
-!
-END DO
-!
-!---------------------------------------------------------------------------
-!
-!*          6. SAVES THE REDELSPERGER NUMBERS
-!              ------------------------------
-!
-IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
-  !
-  ! stores the RED_TH1
-  TZFIELD%CMNHNAME   = 'RED_TH1'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'RED_TH1'
-  TZFIELD%CUNITS     = '1'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_RED_TH1'
-  TZFIELD%NGRID      = 4
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PREDTH1)
-  !
-  ! stores the RED_R1
-  TZFIELD%CMNHNAME   = 'RED_R1'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'RED_R1'
-  TZFIELD%CUNITS     = '1'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_RED_R1'
-  TZFIELD%NGRID      = 4
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PREDR1)
-  !
-  ! stores the RED2_TH3
-  TZFIELD%CMNHNAME   = 'RED2_TH3'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'RED2_TH3'
-  TZFIELD%CUNITS     = '1'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_RED2_TH3'
-  TZFIELD%NGRID      = 4
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PRED2TH3)
-  !
-  ! stores the RED2_R3
-  TZFIELD%CMNHNAME   = 'RED2_R3'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'RED2_R3'
-  TZFIELD%CUNITS     = '1'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_RED2_R3'
-  TZFIELD%NGRID      = 4
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PRED2R3)
-  !
-  ! stores the RED2_THR3
-  TZFIELD%CMNHNAME   = 'RED2_THR3'
-  TZFIELD%CSTDNAME   = ''
-  TZFIELD%CLONGNAME  = 'RED2_THR3'
-  TZFIELD%CUNITS     = '1'
-  TZFIELD%CDIR       = 'XY'
-  TZFIELD%CCOMMENT   = 'X_Y_Z_RED2_THR3'
-  TZFIELD%NGRID      = 4
-  TZFIELD%NTYPE      = TYPEREAL
-  TZFIELD%NDIMS      = 3
-  TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,PRED2THR3)
-  !
-END IF
-!
-!---------------------------------------------------------------------------
-ENDIF ! (Done only if LHARAT is FALSE)
-!
-IF (LHOOK) CALL DR_HOOK('PRANDTL',1,ZHOOK_HANDLE)
-END SUBROUTINE PRANDTL
-END MODULE MODE_PRANDTL
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 9bfa00590f540579dc15c0df3400d9867f71ebc6..e13203dacc08a1ffd063efd0727e0fa47b7ef800 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -5,18 +5,18 @@
 !-----------------------------------------------------------------
       SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,      &
               & KSPLIT,KMODEL_CL,                                     &
-              & OCLOSE_OUT,OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
+              & OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01,    &
               & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HINST_SFU,         &
-              & HMF_UPDRAFT,PIMPL,PTSTEP_UVW, PTSTEP_MET,PTSTEP_SV,   &
-              & HFMFILE,HLUOUT,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,          &
+              & HMF_UPDRAFT,PIMPL,   &
+              & PTSTEP,TPFILE, PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,           &
               & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
               & PRHODJ,PTHVREF,PRHODREF,                              &
               & PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
-              & PPABSM,PUM,PVM,PWM,PTKEM,PSVM,PSRCM,                  &
+              & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
               & PLENGTHM,PLENGTHH,MFMOIST,                            &
               & PBL_DEPTH,PSBL_DEPTH,                                 &
-              & PUT,PVT,PWT,PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,    &
-              & PTHLM,PRM,                                            &
+              & PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT,    &
+              & PTHLT,PRT,                                            &
               & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,              &
               & PHGRAD, PSIGS,                                        &
               & PDRUS_TURB,PDRVS_TURB,                                &
@@ -24,7 +24,7 @@
               & PFLXZTHVMF,PWTH,PWRC,PWSV,PDP,PTP,PTPMF,PTDIFF,PTDISS,&
               & YDDDH,YDLDDH,YDMDDH,                                  &
               & TBUDGETS, KBUDGETS,                                   &
-              & PTR,PDISS,PEDR                                        )
+              & PTR,PDISS,PEDR,PLEM                            )
 !     #################################################################
 !
 !
@@ -222,17 +222,23 @@
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !
-USE MODD_PARAMETERS
+USE MODD_PARAMETERS, ONLY: JPVEXT_TURB
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_CONF
-USE MODD_BUDGET
+USE MODD_BUDGET, ONLY: LBUDGET_U,  LBUDGET_V,  LBUDGET_W,  LBUDGET_TH, LBUDGET_RV, LBUDGET_RC,  &
+                            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, &
+                            TBUDGETDATA
+USE MODD_DYN_n, ONLY : LOCEAN
 USE MODD_FIELD, ONLY: TFIELDDATA,TYPEREAL
 USE MODD_IO, ONLY: TFILEDATA
 USE MODD_LES
+USE MODD_TURB_n, ONLY: XCADAP
 USE MODD_NSV
 !
-USE MODI_BL89
+USE MODE_BL89, ONLY: BL89
 USE MODE_TURB_VER, ONLY : TURB_VER
 !!MODIF AROME
 !USE MODI_ROTATE_WIND
@@ -240,19 +246,22 @@ USE MODE_TURB_VER, ONLY : TURB_VER
 USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
 USE MODI_SHUMAN, ONLY : MZF, MXF, MYF
 USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
 USE MODI_BUDGET_DDH
 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
 !
+USE MODE_BUDGET,         ONLY: BUDGET_STORE_INIT, BUDGET_STORE_END
 USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_SBL
-USE MODE_FMWRIT
+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 DDH_MIX, ONLY  : TYP_DDH
 USE YOMLDDH, ONLY  : TLDDH
@@ -275,8 +284,6 @@ 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)   ::  OCLOSE_OUT   ! switch for syncronous
-                                                      ! file opening
 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
@@ -284,24 +291,18 @@ LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
 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)   ::  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
 CHARACTER(LEN=1),       INTENT(IN)   ::  HINST_SFU    ! temporal location of the
                                                       ! surface friction flux
 REAL,                   INTENT(IN)   ::  PIMPL        ! degree of implicitness
-REAL,                   INTENT(IN)   ::  PTSTEP_UVW   ! Dynamical timestep 
-REAL,                   INTENT(IN)   ::  PTSTEP_MET   ! Timestep for meteorological variables                        
-REAL,                   INTENT(IN)   ::  PTSTEP_SV    ! Timestep for tracer variables
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+REAL,                   INTENT(IN)   ::  PTSTEP       ! timestep 
 !
 CHARACTER(LEN=4),       INTENT(IN)   ::  HMF_UPDRAFT  ! Type of Mass Flux Scheme
-
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDXX,PDYY,PDZZ,PDZX,PDZY
                                         ! metric coefficients
@@ -324,21 +325,20 @@ REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODREF  ! dry density of the
 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 
+! 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) ::  PPABSM      ! Pressure at time t-1
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM,PWM ! wind components
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM        ! passive scal. var.
-REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PSRCM       ! Second-order flux
+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
 !
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PUT,PVT,PWT ! Wind  at t
 !    variables for cloud mixing length
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PCEI ! Cloud Entrainment instability
                                                  ! index to emphasize localy 
@@ -348,9 +348,9 @@ REAL, INTENT(IN)      ::  PCEI_MAX ! maximum threshold for the instability index
 REAL, INTENT(IN)      ::  PCOEF_AMPL_SAT ! saturation of the amplification coefficient
 !
 !   thermodynamical variables which are transformed in conservative var.
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PTHLM       ! conservative pot. temp.
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRM         ! water var.  where 
-                             ! PRM(:,:,:,1) is the conservative mixing ratio        
+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
@@ -361,7 +361,7 @@ 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 
+! saturation
 REAL, DIMENSION(:,:,:,:), INTENT(IN)    ::  PHGRAD
 REAL, DIMENSION(:,:,:), INTENT(OUT)     ::  PSIGS
 REAL, DIMENSION(:,:,:), INTENT(OUT)     ::  PDRUS_TURB   ! evolution of rhoJ*U   by turbulence only
@@ -395,19 +395,20 @@ 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  :: PEDR  ! EDR
+REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL  :: PLEM  ! Mixing length
 !
 !
 !-------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
 !
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::     &
+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
-          ZLM,                        &  ! Turbulent mixing length
+          ZLM,ZLMW,                   &  ! Turbulent mixing length (+ work array)
           ZLEPS,                      &  ! Dissipative length
           ZTRH,                       &  ! 
           ZATHETA,ZAMOIST,            &  ! coefficients for s = f (Thetal,Rnp)
@@ -416,9 +417,9 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::     &
           ZMWTH,ZMWR,ZMTH2,ZMR2,ZMTHR,&  ! 3rd order moments
           ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,&  ! opposite of verticale derivate of 3rd order moments
           ZTHLM                          ! initial potential temp.
-REAL, DIMENSION(SIZE(PRM,1),SIZE(PRM,2),SIZE(PRM,3),SIZE(PRM,4)) ::     &
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),SIZE(PRT,4)) ::     &
           ZRM                            ! initial mixing ratio 
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2)) ::  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,  &
@@ -447,11 +448,9 @@ INTEGER             :: IKTB,IKTE    ! start, end of k loops in physical domain
 INTEGER             :: JRR,JK,JSV   ! loop counters
 INTEGER             :: JI,JJ        ! loop counters
 REAL                :: ZL0          ! Max. Mixing Length in Blakadar formula
-REAL                :: ZALPHA       ! proportionnality constant between Dz/2 and 
-!                                   ! BL89 mixing length near the surface
-!
-!
-TYPE(TFILEDATA) :: TPFILE ! File type to write fields for MesoNH
+REAL                :: ZALPHA       ! work coefficient : 
+                                    ! - proportionnality constant between Dz/2 and 
+!                                   !   BL89 mixing length near the surface
 !
 REAL :: ZTIME1, ZTIME2
 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3))::  ZSHEAR, ZDUDZ, ZDVDZ
@@ -465,6 +464,7 @@ TYPE(TFIELDDATA) :: TZFIELD
 !
 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
@@ -472,7 +472,7 @@ IF (LHARAT .AND. LLES_CALL) THEN
   CALL ABOR1('LHARATU not implemented for option LLES_CALL')
 ENDIF
 
-IKT=SIZE(PTHLM,3)          
+IKT=SIZE(PTHLT,3)          
 IKTB=1+JPVEXT_TURB              
 IKTE=IKT-JPVEXT_TURB
 IKB=KKA+JPVEXT_TURB*KKL
@@ -482,8 +482,11 @@ ZEXPL = 1.- PIMPL
 ZRVORD= XRV / XRD
 !
 !
-ZTHLM(:,:,:) = PTHLM(:,:,:)
-ZRM(:,:,:,:) = PRM(:,:,:,:)
+!Copy data into ZTHLM and ZRM only if needed
+IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. ORMC01) THEN
+  ZTHLM(:,:,:) = PTHLT(:,:,:)
+  ZRM(:,:,:,:) = PRT(:,:,:,:)
+END IF
 !
 !
 !
@@ -494,20 +497,24 @@ ZRM(:,:,:,:) = PRM(:,:,:,:)
 !
 !*      2.1 Cph at t
 !
-ZCP=XCPD
+ZCP(:,:,:)=XCPD
 !
-IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PRM(:,:,:,1)
+IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PRT(:,:,:,1)
 DO JRR = 2,1+KRRL                          ! loop on the liquid components  
-  ZCP(:,:,:)  = ZCP(:,:,:) + XCL * PRM(:,:,:,JRR)
+  ZCP(:,:,:)  = ZCP(:,:,:) + XCL * PRT(:,:,:,JRR)
 END DO
 !
 DO JRR = 2+KRRL,1+KRRL+KRRI                ! loop on the solid components   
-  ZCP(:,:,:)  = ZCP(:,:,:)  + XCI * PRM(:,:,:,JRR)
+  ZCP(:,:,:)  = ZCP(:,:,:)  + XCI * PRT(:,:,:,JRR)
 END DO
 !
 !*      2.2 Exner function at t
 !
-ZEXN(:,:,:) = (PPABSM(:,:,:)/XP00) ** (XRD/XCPD)
+IF (LOCEAN) THEN
+  ZEXN(:,:,:) = 1.
+ELSE
+  ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD)
+END IF
 !
 !*      2.3 dissipative heating coeff a t
 !
@@ -522,23 +529,23 @@ IF (KRRL >=1) THEN
 !
 !*      2.4 Temperature at t
 !
-  ZT(:,:,:) =  PTHLM(:,:,:) * ZEXN(:,:,:)
+  ZT(:,:,:) =  PTHLT(:,:,:) * ZEXN(:,:,:)
 !
 !*       2.5 Lv/Cph/Exn
 !
   IF ( KRRI >= 1 ) THEN 
-    ALLOCATE(ZLVOCPEXNM(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
-    ALLOCATE(ZLSOCPEXNM(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
-    ALLOCATE(ZAMOIST_ICE(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
-    ALLOCATE(ZATHETA_ICE(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
+    ALLOCATE(ZLVOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
+    ALLOCATE(ZLSOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
+    ALLOCATE(ZAMOIST_ICE(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
+    ALLOCATE(ZATHETA_ICE(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
 !
     CALL COMPUTE_FUNCTION_THERMO(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
                                  ZLVOCPEXNM,ZAMOIST,ZATHETA)
     CALL COMPUTE_FUNCTION_THERMO(XALPI,XBETAI,XGAMI,XLSTT,XCI,ZT,ZEXN,ZCP, &
                                  ZLSOCPEXNM,ZAMOIST_ICE,ZATHETA_ICE)
 !
-    WHERE(PRM(:,:,:,2)+PRM(:,:,:,4)>0.0)
-      ZFRAC_ICE(:,:,:) = PRM(:,:,:,4) / ( PRM(:,:,:,2)+PRM(:,:,:,4) )
+    WHERE(PRT(:,:,:,2)+PRT(:,:,:,4)>0.0)
+      ZFRAC_ICE(:,:,:) = PRT(:,:,:,4) / ( PRT(:,:,:,2)+PRT(:,:,:,4) )
     END WHERE
 !
     ZLOCPEXNM(:,:,:) = (1.0-ZFRAC_ICE(:,:,:))*ZLVOCPEXNM(:,:,:) &
@@ -567,7 +574,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   = ''
@@ -579,7 +586,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
@@ -589,21 +596,21 @@ END IF              ! loop end on KRRL >= 1
 ! computes conservative variables
 !
 IF ( KRRL >= 1 ) THEN
-  IF ( KRRI >= 1 ) THEN 
-    ! Rnp at t-1
-    PRM(:,:,:,1)  = PRM(:,:,:,1)  + PRM(:,:,:,2)  + PRM(:,:,:,4)
+  IF ( KRRI >= 1 ) THEN
+    ! Rnp at t
+    PRT(:,:,:,1)  = PRT(:,:,:,1)  + PRT(:,:,:,2)  + PRT(:,:,:,4)
     PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) + PRRS(:,:,:,4)
-    ! Theta_l at t-1
-    PTHLM(:,:,:)  = PTHLM(:,:,:)  - ZLVOCPEXNM(:,:,:) * PRM(:,:,:,2) &
-                                  - ZLSOCPEXNM(:,:,:) * PRM(:,:,:,4)
+    ! Theta_l at t
+    PTHLT(:,:,:)  = PTHLT(:,:,:)  - ZLVOCPEXNM(:,:,:) * PRT(:,:,:,2) &
+                                  - ZLSOCPEXNM(:,:,:) * PRT(:,:,:,4)
     PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZLVOCPEXNM(:,:,:) * PRRS(:,:,:,2) &
                                   - ZLSOCPEXNM(:,:,:) * PRRS(:,:,:,4)
   ELSE
-    ! Rnp at t-1
-    PRM(:,:,:,1)  = PRM(:,:,:,1)  + PRM(:,:,:,2) 
+    ! Rnp at t
+    PRT(:,:,:,1)  = PRT(:,:,:,1)  + PRT(:,:,:,2) 
     PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2)
-    ! Theta_l at t-1
-    PTHLM(:,:,:)  = PTHLM(:,:,:)  - ZLOCPEXNM(:,:,:) * PRM(:,:,:,2)
+    ! Theta_l at t
+    PTHLT(:,:,:)  = PTHLT(:,:,:)  - ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
     PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2)
   END IF
 END IF
@@ -629,21 +636,47 @@ SELECT CASE (HTURBLEN)
 
   CASE ('BL89')
     ZSHEAR=0.
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZSHEAR,ZLM)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM)
 !
-!*      3.2 Delta mixing length
+!*      3.2 RM17 mixing length
+!           ------------------
+
+  CASE ('RM17')
+    ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ,KKA,KKU,KKL),KKA,KKU,KKL))
+    ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ,KKA,KKU,KKL),KKA,KKU,KKL))
+    ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM)
+!
+!*      3.3 Grey-zone combined RM17 & Deardorff mixing lengths 
+!           --------------------------------------------------
+
+  CASE ('ADAP')
+    ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ,KKA,KKU,KKL),KKA,KKU,KKL))
+    ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ,KKA,KKU,KKL),KKA,KKU,KKL))
+    ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+    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.
+    ! For large horizontal grid meshes, this is equal to RM17
+    ! 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.
+    ZLM = MIN(ZLM,XCADAP*ZLMW)
+!
+!*      3.4 Delta mixing length
 !           -------------------
 !
   CASE ('DELT')
-    CALL DELT(ZLM)
+    CALL DELT(PLEM,ODZ=.TRUE.)
 !
-!*      3.3 Deardorff mixing length
+!*      3.5 Deardorff mixing length
 !           -----------------------
 !
   CASE ('DEAR')
     CALL DEAR(ZLM)
 !
-!*      3.4 Blackadar mixing length
+!*      3.6 Blackadar mixing length
 !           -----------------------
 !
   CASE ('BLKR')
@@ -675,26 +708,29 @@ ENDIF  ! end LHARRAT
 !           ------------------
 
 IF (LHARAT) THEN
-ZLEPS=PLENGTHM*(3.75**2.)
+  ZLEPS=PLENGTHM*(3.75**2.)
 ELSE
-ZLEPS=ZLM
+  ZLEPS=ZLM
 ENDIF
 !
 !*      3.7 Correction in the Surface Boundary Layer (Redelsperger 2001)
 !           ----------------------------------------
 !
 ZLMO=XUNDEF
- IF (ORMC01) THEN
-   ZUSTAR=(PSFU**2+PSFV**2)**(0.25)
-    IF (KRR>0) THEN
-     ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV)
-    ELSE
-     ZRVM=0.
-     ZSFRV=0.
-     ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRVM,PSFTH,ZSFRV)
-    END IF
+IF (ORMC01) THEN
+  ZUSTAR=(PSFU**2+PSFV**2)**(0.25)
+  IF (KRR>0) THEN
+    ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV)
+  ELSE
+    ZRVM=0.
+    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,ZLM,ZLEPS)
- END IF
+END IF
+!
+!RMC01 is only applied on RM17 in ADAP
+IF (HTURBLEN=='ADAP') ZLEPS = MIN(ZLEPS,ZLMW*XCADAP)
 !
 !*      3.8 Mixing length in external points (used if HTURBDIM="3DIM")
 !           ----------------------------------------------------------
@@ -715,8 +751,8 @@ IF ( HINST_SFU == 'T' ) THEN
 !
 !
   IF (CPROGRAM=='AROME ') THEN
-    ZUSLOPE=PUM(:,:,KKA)
-    ZVSLOPE=PVM(:,:,KKA)
+    ZUSLOPE=PUT(:,:,KKA)
+    ZVSLOPE=PVT(:,:,KKA)
   ELSE
 !    CALL ROTATE_WIND(PUT,PVT,PWT,                       &
 !                     PDIRCOSXW, PDIRCOSYW, PDIRCOSZW,   &
@@ -737,7 +773,7 @@ IF ( HINST_SFU == 'T' ) THEN
 !*      4.3 rotate the wind at time t-delta t
 !
   IF (CPROGRAM/='AROME ') THEN
-!    CALL ROTATE_WIND(PUM,PVM,PWM,                       &
+!    CALL ROTATE_WIND(PUT,PVT,PWT,                       &
 !                     PDIRCOSXW, PDIRCOSYW, PDIRCOSZW,   &
 !                     PCOSSLOPE,PSINSLOPE,               &
 !                     PDXX,PDYY,PDZZ,                    &
@@ -751,11 +787,11 @@ ELSE
 !*      4.4 rotate the wind at time t-delta t
 !
   IF (CPROGRAM=='AROME ') THEN
-    ZUSLOPE=PUM(:,:,KKA)
-    ZVSLOPE=PVM(:,:,KKA)
+    ZUSLOPE=PUT(:,:,KKA)
+    ZVSLOPE=PVT(:,:,KKA)
   ELSE
 !
-!    CALL ROTATE_WIND(PUM,PVM,PWM,                       &
+!    CALL ROTATE_WIND(PUT,PVT,PWT,                       &
 !                     PDIRCOSXW, PDIRCOSYW, PDIRCOSZW,   &
 !                     PCOSSLOPE,PSINSLOPE,               &
 !                     PDXX,PDYY,PDZZ,                    &
@@ -775,8 +811,8 @@ END IF
 !
 ZTAU11M(:,:) =2./3.*(  (1.+ (PZZ (:,:,IKB+KKL)-PZZ (:,:,IKB))  &
                            /(PDZZ(:,:,IKB+KKL)+PDZZ(:,:,IKB))  &
-                       )   *PTKEM(:,:,IKB)                   &
-                     -0.5  *PTKEM(:,:,IKB+KKL)                 &
+                       )   *PTKET(:,:,IKB)                   &
+                     -0.5  *PTKET(:,:,IKB+KKL)                 &
                     )
 ZTAU12M(:,:) =0.0
 ZTAU22M(:,:) =ZTAU11M(:,:)
@@ -792,129 +828,217 @@ ZMTH2 = 0.     ! w'th'2
 ZMR2  = 0.     ! w'r'2
 ZMTHR = 0.     ! w'th'r'
 
-IF (HTOM=='TM06') CALL TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
-!
-ZFWTH = -GZ_M_W(ZMWTH,PDZZ, KKA, KKU, KKL)    ! -d(w'2th' )/dz
-ZFWR  = -GZ_M_W(ZMWR, PDZZ, KKA, KKU, KKL)    ! -d(w'2r'  )/dz
-ZFTH2 = -GZ_W_M(ZMTH2,PDZZ, KKA, KKU, KKL)    ! -d(w'th'2 )/dz
-ZFR2  = -GZ_W_M(ZMR2, PDZZ, KKA, KKU, KKL)    ! -d(w'r'2  )/dz
-ZFTHR = -GZ_W_M(ZMTHR,PDZZ, KKA, KKU, KKL)    ! -d(w'th'r')/dz
-!
-ZFWTH(:,:,IKTE:) = 0.
-ZFWTH(:,:,:IKTB) = 0.
-ZFWR (:,:,IKTE:) = 0.
-ZFWR (:,:,:IKTB) = 0.
-ZFTH2(:,:,IKTE:) = 0.
-ZFTH2(:,:,:IKTB) = 0.
-ZFR2 (:,:,IKTE:) = 0.
-ZFR2 (:,:,:IKTB) = 0.
-ZFTHR(:,:,IKTE:) = 0.
-ZFTHR(:,:,:IKTB) = 0.
+IF (HTOM=='TM06') THEN
+  CALL TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
+!
+  ZFWTH = -GZ_M_W(ZMWTH,PDZZ, KKA,KKU,KKL)    ! -d(w'2th' )/dz
+  !ZFWR  = -GZ_M_W(ZMWR, PDZZ,KKA,KKU,KKL)    ! -d(w'2r'  )/dz
+  ZFTH2 = -GZ_W_M(ZMTH2,PDZZ,KKA,KKU,KKL)    ! -d(w'th'2 )/dz
+  !ZFR2  = -GZ_W_M(ZMR2, PDZZ,KKA,KKU,KKL)    ! -d(w'r'2  )/dz
+  !ZFTHR = -GZ_W_M(ZMTHR,PDZZ,KKA,KKU,KKL)    ! -d(w'th'r')/dz
+!
+  ZFWTH(:,:,IKTE:) = 0.
+  ZFWTH(:,:,:IKTB) = 0.
+  !ZFWR (:,:,IKTE:) = 0.
+  !ZFWR (:,:,:IKTB) = 0.
+  ZFWR  = 0.
+  ZFTH2(:,:,IKTE:) = 0.
+  ZFTH2(:,:,:IKTB) = 0.
+  !ZFR2 (:,:,IKTE:) = 0.
+  !ZFR2 (:,:,:IKTB) = 0.
+  ZFR2  = 0.
+  !ZFTHR(:,:,IKTE:) = 0.
+  !ZFTHR(:,:,:IKTB) = 0.
+  ZFTHR = 0.
+ELSE
+  ZFWTH = 0.
+  ZFWR  = 0.
+  ZFTH2 = 0.
+  ZFR2  = 0.
+  ZFTHR = 0.
+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_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_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
+
 CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI,               &
           OTURB_FLX,                                     &
           HTURBDIM,HTOM,PIMPL,ZEXPL,                     &
-          PTSTEP_MET,TPFILE,                                 &
+          PTSTEP,TPFILE,                                 &
           PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,        &
           PCOSSLOPE,PSINSLOPE,                           &
           PRHODJ,PTHVREF,                                &
           PSFTH,PSFRV,PSFSV,PSFTH,PSFRV,PSFSV,           &
           ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU33M,               &
-          PUM,PVM,PWM,ZUSLOPE,ZVSLOPE,PTHLM,PRM,PSVM,    &
-          PTKEM,ZLM,PLENGTHM,PLENGTHH,ZLEPS,MFMOIST,     &
-          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCM,ZFRAC_ICE,     &
+          PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,    &
+          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,              &
           PDP,PTP,PSIGS,PWTH,PWRC,PWSV                   )
-!
 
-IF (LBUDGET_U) CALL BUDGET_DDH (PRUS,1,'VTURB_BU_RU',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_V) CALL BUDGET_DDH (PRVS,2,'VTURB_BU_RV',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_W) CALL BUDGET_DDH (PRWS,3,'VTURB_BU_RW',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_TH)  THEN
-  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) + ZLSOCPEXNM * PRRS(:,:,:,4),4,'VTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
-  ELSE IF ( KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2),4,'VTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+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_DDH (PRTHLS,4,'VTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:, :, :) )
   END IF
 END IF
-IF (LBUDGET_SV) THEN
-  DO JSV = 1,NSV
-    CALL BUDGET_DDH (PRSVS(:,:,:,JSV),JSV+12,'VTURB_BU_RSV',YDDDH, YDLDDH, YDMDDH)
+
+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_SV )  THEN
+  DO JSV = 1, NSV
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:, :, :, JSV) )
   END DO
 END IF
-IF (LBUDGET_RV) THEN
-  IF ( KRRI >= 1 .AND. KRRL >= 1) THEN
-    CALL BUDGET_DDH (PRRS(:,:,:,1)-PRRS(:,:,:,2)-PRRS(:,:,:,4),6,'VTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
-  ELSE IF ( KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRRS(:,:,:,1)-PRRS(:,:,:,2),6,'VTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
-  ELSE 
-    CALL BUDGET_DDH (PRRS(:,:,:,1),6,'VTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
-  END IF
-END IF  
-IF (LBUDGET_RC) CALL BUDGET_DDH (PRRS(:,:,:,2),7,'VTURB_BU_RRC',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_RI) CALL BUDGET_DDH (PRRS(:,:,:,4),9,'VTURB_BU_RRI',YDDDH, YDLDDH, YDMDDH)
-!
 !
-IF (HTURBDIM=='3DIM') THEN
-!!!!MODIF AROME
-!  CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP_UVW,      &
-!          PTSTEP_MET, PTSTEP_SV, HLBCX,HLBCY,                  &
-!          OCLOSE_OUT,OTURB_FLX,OSUBG_COND,                     &
-!          HFMFILE,HLUOUT,                                      &
+#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_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_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,                    &
+!          TPFILE,                                              &
 !          PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,                        &
 !          PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,                       &
 !          PCOSSLOPE,PSINSLOPE,                                 &
 !          PRHODJ,PTHVREF,                                      &
 !          PSFTH,PSFRV,PSFSV,                                   &
 !          ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M,             &
-!          PUM,PVM,PWM,ZUSLOPE,ZVSLOPE,PTHLM,PRM,PSVM,          &
-!          PTKEM,ZLM,ZLEPS,                                     &
-!          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCM,ZFRAC_ICE,           &
-!          ZDP,ZTP,PSIGS,                                       &
-!          PHGRAD,                                              &
+!          PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,          &
+!          PTKET,PLEM,ZLEPS,                                    &
+!          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,           &
+!          PDYP,PTHP,PSIGS,                                     &
 !          ZTRH,                                                &
 !          PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS                     )
-END IF
-!
-!
-IF (LBUDGET_U) CALL BUDGET_DDH (PRUS,1,'HTURB_BU_RU',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_V) CALL BUDGET_DDH (PRVS,2,'HTURB_BU_RV',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_W) CALL BUDGET_DDH (PRWS,3,'HTURB_BU_RW',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_TH)  THEN
-  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) + ZLSOCPEXNM * PRRS(:,:,:,4),4,'HTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
-  ELSE IF ( KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2),4,'HTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
-  ELSE
-    CALL BUDGET_DDH (PRTHLS,4,'HTURB_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+
+  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
-END IF
-IF (LBUDGET_SV) THEN
-  DO JSV = 1,NSV
-    CALL BUDGET_DDH (PRSVS(:,:,:,JSV),JSV+12,'HTURB_BU_RSV',YDDDH, YDLDDH, YDMDDH)
-  END DO
-END IF
-IF (LBUDGET_RV) THEN
-  IF ( KRRI >= 1 .AND. KRRL >= 1) THEN
-    CALL BUDGET_DDH (PRRS(:,:,:,1)-PRRS(:,:,:,2)-PRRS(:,:,:,4),6,'HTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
-  ELSE IF ( KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRRS(:,:,:,1)-PRRS(:,:,:,2),6,'HTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
-  ELSE 
-    CALL BUDGET_DDH (PRRS(:,:,:,1),6,'HTURB_BU_RRV',YDDDH, YDLDDH, YDMDDH)
+
+  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
-END IF  
-IF (LBUDGET_RC) CALL BUDGET_DDH (PRRS(:,:,:,2),7,'HTURB_BU_RRC',YDDDH, YDLDDH, YDMDDH)
-IF (LBUDGET_RI) CALL BUDGET_DDH (PRRS(:,:,:,4),9,'HTURB_BU_RRI',YDDDH, YDLDDH, YDMDDH)
-!
+
+  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
+#ifdef REPRO48
+#else
+END IF
+#endif
 !----------------------------------------------------------------------------
 !
 !*      6. EVOLUTION OF THE TKE AND ITS DISSIPATION 
@@ -929,10 +1053,21 @@ IF (LBUDGET_RI) CALL BUDGET_DDH (PRRS(:,:,:,4),9,'HTURB_BU_RRI',YDDDH, YDLDDH, Y
 !  6.2 TKE evolution equation
 
 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,PTKEM,ZLM,ZLEPS,PDP,ZTRH,       &
+CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,       &
                    & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,            &
-                   & PTSTEP_MET,PIMPL,ZEXPL,                         &
+                   & PTSTEP,PIMPL,ZEXPL,                         &
                    & HTURBLEN,HTURBDIM,                              &
                    & TPFILE,OTURB_DIAG,           &
                    & PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,&
@@ -940,14 +1075,15 @@ CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,ZLM,ZLEPS,PDP,ZTRH,       &
                    & PEDR=PEDR)
 IF (LBUDGET_TH)  THEN
   IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) + ZLSOCPEXNM * PRRS(:,:,:,4),4,'DISSH_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
+                                                          & + ZLSOCPEXNM * PRRS(:,:,:,4) )
   ELSE IF ( KRRL >= 1 ) THEN
-    CALL BUDGET_DDH (PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2),4,'DISSH_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
   ELSE
-    CALL BUDGET_DDH (PRTHLS,4,'DISSH_BU_RTH',YDDDH, YDLDDH, YDMDDH)
+    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
   END IF
 END IF
-
+!
 ENDIF
 !
 !----------------------------------------------------------------------------
@@ -969,7 +1105,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
   TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,ZLM)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
 !
   IF (KRR /= 0) THEN
 !
@@ -985,7 +1121,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,PTHLM)
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PTHLT)
 !
 ! stores the conservative mixing ratio
 !
@@ -999,7 +1135,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 3
     TZFIELD%LTIMEDEP   = .TRUE.
-    CALL IO_Field_write(TPFILE,TZFIELD,PRM(:,:,:,1))
+    CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PRT(:,:,:,1))
    END IF
 END IF
 !
@@ -1016,23 +1152,23 @@ PDRSVS_TURB  = PRSVS - PDRSVS_TURB
 !
 IF ( KRRL >= 1 ) THEN
   IF ( KRRI >= 1 ) THEN
-    PRM(:,:,:,1)  = PRM(:,:,:,1)  - PRM(:,:,:,2)  - PRM(:,:,:,4)
+    PRT(:,:,:,1)  = PRT(:,:,:,1)  - PRT(:,:,:,2)  - PRT(:,:,:,4)
     PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2) - PRRS(:,:,:,4)
-    PTHLM(:,:,:)  = PTHLM(:,:,:)  + ZLVOCPEXNM(:,:,:) * PRM(:,:,:,2) &
-                                  + ZLSOCPEXNM(:,:,:) * PRM(:,:,:,4)
+    PTHLT(:,:,:)  = PTHLT(:,:,:)  + ZLVOCPEXNM(:,:,:) * PRT(:,:,:,2) &
+                                  + ZLSOCPEXNM(:,:,:) * PRT(:,:,:,4)
     PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLVOCPEXNM(:,:,:) * PRRS(:,:,:,2) &
                                   + ZLSOCPEXNM(:,:,:) * PRRS(:,:,:,4)
 !
     DEALLOCATE(ZLVOCPEXNM)
     DEALLOCATE(ZLSOCPEXNM)
   ELSE
-    PRM(:,:,:,1)  = PRM(:,:,:,1)  - PRM(:,:,:,2) 
+    PRT(:,:,:,1)  = PRT(:,:,:,1)  - PRT(:,:,:,2) 
     PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2)
-    PTHLM(:,:,:)  = PTHLM(:,:,:)  + ZLOCPEXNM(:,:,:) * PRM(:,:,:,2)
+    PTHLT(:,:,:)  = PTHLT(:,:,:)  + ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
     PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2)
   END IF
 END IF
-!
+
 !----------------------------------------------------------------------------
 !
 !*      9. LES averaged surface fluxes
@@ -1067,16 +1203,16 @@ IF (LLES_CALL) THEN
 !          ------------------------------------------------
 !
   IF (HTURBDIM=="1DIM") THEN
-    CALL LES_MEAN_SUBGRID(2./3.*PTKEM,X_LES_SUBGRID_U2)
-    CALL LES_MEAN_SUBGRID(2./3.*PTKEM,X_LES_SUBGRID_V2)
-    CALL LES_MEAN_SUBGRID(2./3.*PTKEM,X_LES_SUBGRID_W2)
-    CALL LES_MEAN_SUBGRID(2./3.*PTKEM*MZF(GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL),&
+    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_U2)
+    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_V2)
+    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_W2)
+    CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(PTHLT,PDZZ, KKA, KKU, KKL),&
                           KKA, KKU, KKL),X_LES_RES_ddz_Thl_SBG_W2)
     IF (KRR>=1) &
-    CALL LES_MEAN_SUBGRID(2./3.*PTKEM*MZF(GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL),&
+    CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(PRT(:,:,:,1),PDZZ, KKA, KKU, KKL),&
                          &KKA, KKU, KKL),X_LES_RES_ddz_Rt_SBG_W2)
     DO JSV=1,NSV
-      CALL LES_MEAN_SUBGRID(2./3.*PTKEM*MZF(GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL), &
+      CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(PSVT(:,:,:,JSV),PDZZ, KKA, KKU, KKL), &
                            &KKA, KKU, KKL), X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
     END DO
   END IF
@@ -1091,13 +1227,14 @@ IF (LLES_CALL) THEN
 !
 !* presso-correlations for subgrid Tke are equal to zero.
 !
-  ZLM = 0.
-  CALL LES_MEAN_SUBGRID(ZLM,X_LES_SUBGRID_WP)
+  ZLEPS = 0. !ZLEPS is used as a work array (not used anymore)
+  CALL LES_MEAN_SUBGRID(ZLEPS,X_LES_SUBGRID_WP)
 !
   CALL SECOND_MNH(ZTIME2)
   XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 END IF
 !
+IF(PRESENT(PLEM)) PLEM = ZLM
 !----------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
@@ -1209,7 +1346,7 @@ IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments 
 !
-REAL                                  :: PALP,PBETA,PGAM,PLTT,PC
+REAL,                   INTENT(IN)    :: PALP,PBETA,PGAM,PLTT,PC
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PT,PEXN,PCP
 !
 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLOCPEXN
@@ -1237,7 +1374,7 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
 !
 !*      1.3 saturation  mixing ratio at t
 !
-  ZRVSAT(:,:,:) =  ZRVSAT(:,:,:) * ZEPS / ( PPABSM(:,:,:) - ZRVSAT(:,:,:) )
+  ZRVSAT(:,:,:) =  ZRVSAT(:,:,:) * ZEPS / ( PPABST(:,:,:) - ZRVSAT(:,:,:) )
 !
 !*      1.4 compute the saturation mixing ratio derivative (rvs')
 !
@@ -1251,7 +1388,7 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
 !*      1.6 compute Atheta
 !
   PATHETA(:,:,:)= PAMOIST(:,:,:) * PEXN(:,:,:) *                             &
-        ( ( ZRVSAT(:,:,:) - PRM(:,:,:,1) ) * PLOCPEXN(:,:,:) /               &
+        ( ( ZRVSAT(:,:,:) - PRT(:,:,:,1) ) * PLOCPEXN(:,:,:) /               &
           ( 1. + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) )        *               &
           (                                                                  &
            ZRVSAT(:,:,:) * (1. + ZRVSAT(:,:,:)/ZEPS)                         &
@@ -1270,7 +1407,7 @@ IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE)
 END SUBROUTINE COMPUTE_FUNCTION_THERMO
 !
 !     ####################
-      SUBROUTINE DELT(PLM)
+      SUBROUTINE DELT(PLM,ODZ)
 !     ####################
 !!
 !!****  *DELT* routine to compute mixing length for DELT case
@@ -1292,6 +1429,7 @@ END SUBROUTINE COMPUTE_FUNCTION_THERMO
 !*       0.1   Declarations of dummy arguments 
 !
 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLM
+LOGICAL,                INTENT(IN)    :: ODZ
 !
 !*       0.2   Declarations of local variables
 !
@@ -1301,16 +1439,29 @@ REAL                :: ZD           ! distance to the surface
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE)
-DO JK = IKTB,IKTE ! 1D turbulence scheme
-  PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK)
-END DO
-PLM(:,:,KKU) = PLM(:,:,IKE)
-PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
-IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
-  IF ( L2D) THEN
-    PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) 
-  ELSE
-    PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
+IF (ODZ) THEN
+  ! Dz is take into account in the computation
+  DO JK = IKTB,IKTE ! 1D turbulence scheme
+    PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK)
+  END DO
+  PLM(:,:,KKU) = PLM(:,:,IKE)
+  PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
+  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
+    IF ( L2D) THEN
+      PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) 
+    ELSE
+      PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
+    END IF
+  END IF
+ELSE
+  ! Dz not taken into account in computation to assure invariability with vertical grid mesh
+  PLM=1.E10
+  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
+    IF ( L2D) THEN
+      PLM(:,:,:) = MXF(PDXX(:,:,:))
+    ELSE
+      PLM(:,:,:) = (MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./2.)
+    END IF
   END IF
 END IF
 !
@@ -1320,17 +1471,28 @@ END IF
 IF (.NOT. ORMC01) THEN
   ZALPHA=0.5**(-1.5)
   !
-  DO JJ=1,SIZE(PUM,2)
-    DO JI=1,SIZE(PUM,1)
-      DO JK=IKTB,IKTE
-        ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
-        -PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
-        IF ( PLM(JI,JJ,JK)>ZD) THEN
-          PLM(JI,JJ,JK)=ZD
-        ELSE
-          EXIT
-        ENDIF
-      END DO
+  DO JJ=1,SIZE(PUT,2)
+    DO JI=1,SIZE(PUT,1)
+      IF (LOCEAN) THEN
+        DO JK=IKTE,IKTB,-1
+          ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
+          IF ( PLM(JI,JJ,JK)>ZD) THEN
+            PLM(JI,JJ,JK)=ZD
+          ELSE
+            EXIT
+          ENDIF
+       END DO
+      ELSE
+        DO JK=IKTB,IKTE
+          ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
+          -PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
+          IF ( PLM(JI,JJ,JK)>ZD) THEN
+            PLM(JI,JJ,JK)=ZD
+          ELSE
+            EXIT
+          ENDIF
+        END DO
+      ENDIF   
     END DO
   END DO
 END IF
@@ -1345,7 +1507,7 @@ END SUBROUTINE DELT
       SUBROUTINE DEAR(PLM)
 !     ####################
 !!
-!!****  *DELT* routine to compute mixing length for DEARdorff case
+!!****  *DEAR* routine to compute mixing length for DEARdorff case
 !
 !!    AUTHOR
 !!    ------
@@ -1370,9 +1532,10 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLM
 !*       0.2   Declarations of local variables
 !
 REAL                :: ZD           ! distance to the surface
-REAL, DIMENSION(:,:), ALLOCATABLE  ::   ZWORK2D
+REAL                :: ZVAR         ! Intermediary variable
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) ::   ZWORK2D
 !
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::     &
+REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) ::     &
             ZDTHLDZ,ZDRTDZ,     &!dtheta_l/dz, drt_dz used for computing the stablity
 !                                ! criterion 
             ZETHETA,ZEMOIST             !coef ETHETA and EMOIST
@@ -1383,9 +1546,8 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::     &
 !   initialize the mixing length with the mesh grid
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE)
-DO JK = IKTB,IKTE ! 1D turbulence scheme
-  PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK)
-END DO
+! 1D turbulence scheme
+PLM(:,:,IKTB:IKTE) = PZZ(:,:,IKTB+KKL:IKTE+KKL) - PZZ(:,:,IKTB:IKTE)
 PLM(:,:,KKU) = PLM(:,:,IKE)
 PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
 IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
@@ -1397,53 +1559,98 @@ IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
 END IF
 !   compute a mixing length limited by the stability
 !
-ALLOCATE(ZWORK2D(SIZE(PUM,1),SIZE(PUM,2)))
-!
-ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,ZLOCPEXNM,ZATHETA,PSRCM)
-ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,ZLOCPEXNM,ZAMOIST,PSRCM)
+ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
+ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
+!
+IF (KRR>0) THEN
+  DO JK = IKTB+1,IKTE-1
+    DO JJ=1,SIZE(PUT,2)
+      DO JI=1,SIZE(PUT,1)
+        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
+        ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK    ,1))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PRT(JI,JJ,JK    ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK    ))
+        IF (LOCEAN) THEN
+          ZVAR=XG*(XALPHAOC*ZDTHLDZ(JI,JJ,JK)-XBETAOC*ZDRTDZ(JI,JJ,JK))
+        ELSE
+          ZVAR=XG/PTHVREF(JI,JJ,JK)*                                                  &
+             (ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK))
+        END IF
+        !
+        IF (ZVAR>0.) THEN
+          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
+                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
+        END IF
+      END DO
+    END DO
+  END DO
+ELSE! For dry atmos or unsalted ocean runs
+  DO JK = IKTB+1,IKTE-1
+    DO JJ=1,SIZE(PUT,2)
+      DO JI=1,SIZE(PUT,1)
+        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
+        IF (LOCEAN) THEN
+          ZVAR= XG*XALPHAOC*ZDTHLDZ(JI,JJ,JK)
+        ELSE
+          ZVAR= XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
+        END IF
 !
-DO JK = IKTB+1,IKTE-1
-      ZDTHLDZ(:,:,JK)= 0.5*((PTHLM(:,:,JK+KKL)-PTHLM(:,:,JK))/PDZZ(:,:,JK+KKL)+      &
-                         (PTHLM(:,:,JK)-PTHLM(:,:,JK-KKL))/PDZZ(:,:,JK))
-      ZDRTDZ(:,:,JK)= 0.5*((PRM(:,:,JK+KKL,1)-PRM(:,:,JK,1))/PDZZ(:,:,JK+KKL)+      &
-                         (PRM(:,:,JK,1)-PRM(:,:,JK-KKL,1))/PDZZ(:,:,JK))
-       ZWORK2D(:,:)=XG/PTHVREF(:,:,JK)*                                           &
-            (ZETHETA(:,:,JK)*ZDTHLDZ(:,:,JK)+ZEMOIST(:,:,JK)*ZDRTDZ(:,:,JK))
-  !
-  WHERE(ZWORK2D(:,:)>0.)
-    PLM(:,:,JK)=MAX(1.E-10,MIN(PLM(:,:,JK),                &
-                    0.76* SQRT(PTKEM(:,:,JK)/ZWORK2D(:,:))))
-  END WHERE
-END DO
+        IF (ZVAR>0.) THEN
+          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
+                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
+        END IF
+      END DO
+    END DO
+  END DO
+END IF
 !  special case near the surface 
-ZDTHLDZ(:,:,IKB)=(PTHLM(:,:,IKB+KKL)-PTHLM(:,:,IKB))/PDZZ(:,:,IKB+KKL)
-ZDRTDZ(:,:,IKB)=(PRM(:,:,IKB+KKL,1)-PRM(:,:,IKB,1))/PDZZ(:,:,IKB+KKL)
+ZDTHLDZ(:,:,IKB)=(PTHLT(:,:,IKB+KKL)-PTHLT(:,:,IKB))/PDZZ(:,:,IKB+KKL)
+! For dry simulations
+IF (KRR>0) THEN
+  ZDRTDZ(:,:,IKB)=(PRT(:,:,IKB+KKL,1)-PRT(:,:,IKB,1))/PDZZ(:,:,IKB+KKL)
+ELSE
+  ZDRTDZ(:,:,IKB)=0
+ENDIF
 !
-ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)*                                           &
-            (ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
+IF (LOCEAN) THEN
+  ZWORK2D(:,:)=XG*(XALPHAOC*ZDTHLDZ(:,:,IKB)-XBETAOC*ZDRTDZ(:,:,IKB))
+ELSE
+  ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)*                                           &
+              (ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
+END IF
 WHERE(ZWORK2D(:,:)>0.)
-  PLM(:,:,IKB)=MAX(1.E-10,MIN( PLM(:,:,JK),                 &
-                    0.76* SQRT(PTKEM(:,:,IKB)/ZWORK2D(:,:))))
+  PLM(:,:,IKB)=MAX(XMNH_EPSILON,MIN( PLM(:,:,IKB),                 &
+                    0.76* SQRT(PTKET(:,:,IKB)/ZWORK2D(:,:))))
 END WHERE
 !
-DEALLOCATE(ZWORK2D)
-!
 !  mixing length limited by the distance normal to the surface (with the same factor as for BL89)
 !
 IF (.NOT. ORMC01) THEN
   ZALPHA=0.5**(-1.5)
   !
-  DO JJ=1,SIZE(PUM,2)
-    DO JI=1,SIZE(PUM,1)
-      DO JK=IKTB,IKTE
-        ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
-          *PDIRCOSZW(JI,JJ)
-        IF ( PLM(JI,JJ,JK)>ZD) THEN
-          PLM(JI,JJ,JK)=ZD
-        ELSE
-          EXIT
-        ENDIF
-      END DO
+  DO JJ=1,SIZE(PUT,2)
+    DO JI=1,SIZE(PUT,1)
+      IF (LOCEAN) THEN
+        DO JK=IKTE,IKTB,-1
+          ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
+          IF ( PLM(JI,JJ,JK)>ZD) THEN
+            PLM(JI,JJ,JK)=ZD
+          ELSE
+            EXIT
+          ENDIF
+        END DO
+      ELSE
+        DO JK=IKTB,IKTE
+          ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
+            *PDIRCOSZW(JI,JJ)
+          IF ( PLM(JI,JJ,JK)>ZD) THEN
+            PLM(JI,JJ,JK)=ZD
+          ELSE
+            EXIT
+          ENDIF
+        END DO
+      ENDIF 
     END DO
   END DO
 END IF
@@ -1509,10 +1716,10 @@ IMPLICIT NONE
 REAL :: ZPENTE            ! Slope of the amplification straight line
 REAL :: ZCOEF_AMPL_CEI_NUL! Ordonnate at the origin of the
                           ! amplification straight line
-REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: ZCOEF_AMPL
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCOEF_AMPL
                           ! Amplification coefficient of the mixing length
                           ! when the instability criterium is verified 
-REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: ZLM_CLOUD
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
                           ! Turbulent mixing length in the clouds
 !
 !-------------------------------------------------------------------------------
@@ -1552,14 +1759,14 @@ ELSE
 !
 !*         3.1 BL89 mixing length
 !           ------------------
-  CASE ('BL89')
+  CASE ('BL89','RM17','ADAP')
     ZSHEAR=0.
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZSHEAR,ZLM_CLOUD)
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD)
 !
 !*         3.2 Delta mixing length
 !           -------------------
   CASE ('DELT')
-    CALL DELT(ZLM_CLOUD)
+    CALL DELT(ZLM_CLOUD,ODZ=.TRUE.)
 !
 !*         3.3 Deardorff mixing length
 !           -----------------------
@@ -1584,7 +1791,7 @@ IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   TZFIELD%NTYPE      = TYPEREAL
   TZFIELD%NDIMS      = 3
   TZFIELD%LTIMEDEP   = .TRUE.
-  CALL IO_Field_write(TPFILE,TZFIELD,ZLM)
+  CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
 ENDIF
 !
 ! Amplification of the mixing length when the criteria are verified
@@ -1610,7 +1817,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   = ''
@@ -1621,7 +1828,7 @@ 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
 !
diff --git a/src/mesonh/turb/th_r_from_thl_rt_1d.f90 b/src/mesonh/turb/th_r_from_thl_rt_1d.f90
index dbb4d564a66e743470b2204ecd9d0d813b092cc2..8c58fe8d0f9b784c89c6f3f0455e2fe425b6b11d 100644
--- a/src/mesonh/turb/th_r_from_thl_rt_1d.f90
+++ b/src/mesonh/turb/th_r_from_thl_rt_1d.f90
@@ -81,7 +81,7 @@ IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
 !
-CHARACTER(len=1),   INTENT(IN) :: HFRAC_ICE
+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
diff --git a/src/mesonh/turb/th_r_from_thl_rt_2d.f90 b/src/mesonh/turb/th_r_from_thl_rt_2d.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d610a72bda073171f6154bdc8e6a6f1e6763300f
--- /dev/null
+++ b/src/mesonh/turb/th_r_from_thl_rt_2d.f90
@@ -0,0 +1,128 @@
+!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
new file mode 100644
index 0000000000000000000000000000000000000000..0591be50c12cee5be78a7522530ddfa5617a35f1
--- /dev/null
+++ b/src/mesonh/turb/th_r_from_thl_rt_3d.f90
@@ -0,0 +1,126 @@
+!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/turb.f90 b/src/mesonh/turb/turb.f90
index f1b3f9e114dfd662ae1ef65ef1d29d0b8856c700..fff7f9fed9f377f94378cedf8275d8bf8aab960c 100644
--- a/src/mesonh/turb/turb.f90
+++ b/src/mesonh/turb/turb.f90
@@ -239,7 +239,7 @@ END MODULE MODI_TURB
 !!    IMPLICIT ARGUMENTS
 !!    ------------------
 !!
-!!       MODD_PARAMETERS : JPVEXT  number of marginal vertical points
+!!       MODD_PARAMETERS : JPVEXT_TURB  number of marginal vertical points
 !!
 !!       MODD_CONF      : CCONF model configuration (start/restart)
 !!                        L1D   switch for 1D model version
@@ -362,12 +362,11 @@ END MODULE MODI_TURB
 !
 USE PARKIND1, ONLY : JPRB
 USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-!
-use modd_budget,      only: lbudget_u,  lbudget_v,  lbudget_w,  lbudget_th, lbudget_rv, lbudget_rc,  &
-                            lbudget_rr, lbudget_ri, lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv,  &
+USE MODD_BUDGET, ONLY: LBUDGET_U,  LBUDGET_V,  LBUDGET_W,  LBUDGET_TH, LBUDGET_RV, LBUDGET_RC,  &
+                            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, &
-                            TBUDGETDATA
+                            TBUDGETS
 USE MODD_CONF
 USE MODD_CST
 USE MODD_CTURB
@@ -389,7 +388,6 @@ USE MODI_ROTATE_WIND
 USE MODI_TURB_HOR_SPLT 
 USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
 USE MODI_SHUMAN
-USE MODI_GRADIENT_M
 USE MODI_LES_MEAN_SUBGRID
 USE MODI_RMC01
 USE MODI_GRADIENT_W
@@ -397,10 +395,10 @@ USE MODE_TM06, ONLY: TM06
 USE MODI_UPDATE_LM
 USE MODI_GET_HALO
 !
-use mode_budget,         only: Budget_store_init, Budget_store_end
+USE MODE_BUDGET,         ONLY: BUDGET_STORE_INIT, BUDGET_STORE_END
 USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_SBL
-use mode_sources_neg_correct, only: Sources_neg_correct
+USE MODE_SOURES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT
 !
 USE MODI_EMOIST
 USE MODI_ETHETA
@@ -434,11 +432,11 @@ LOGICAL,                INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
 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
+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
+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 
@@ -843,8 +841,6 @@ SELECT CASE (HTURBLEN)
 !
 END SELECT
 !
-!
-!
 !*      3.5 Mixing length modification for cloud
 !           -----------------------
 IF (KMODEL_CL==KMI .AND. HTURBLEN_CL/='NONE') CALL CLOUD_MODIF_LM
@@ -1532,7 +1528,6 @@ ELSE
     END IF
   END IF
 END IF
-
 !
 !  mixing length limited by the distance normal to the surface 
 !  (with the same factor as for BL89)