diff --git a/src/arome/aux/modd_frc.f90 b/src/arome/aux/modd_frc.f90
new file mode 100644
index 0000000000000000000000000000000000000000..e1d77673e8619106b68cb329fa2ed58c96844cf3
--- /dev/null
+++ b/src/arome/aux/modd_frc.f90
@@ -0,0 +1,112 @@
+!MNH_LIC Copyright 1996-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 MODD_FRC
+!     ###############
+!
+!!***  *MODD_FRC -  Declarative module for the forcing fields
+!!
+!!    PURPOSE
+!!    -------
+!       This module contains NFRC 1D-arrays used by FORCING (geostrophic wind
+!     components, large scale vertical wind, theta and humidity profiles when
+!     the relaxation option is used,large scale theta and humidity gradients
+!     and the translation speed of the domain of simulation.
+!     The following control parameters are used by FORCING:
+!     - LGEOST_UV_FRC and LGEOST_TH_FRC
+!     - LTEND_THRV_FRC and LTEND_UV_FRC
+!     - LVERT_MOTION_FRC  
+!     - LRELAX_THRV_FRC, LRELAX_UV_FRC and LRELAX_UVMEAN_FRC using:
+!         XRELAX_TIME_FRC, XRELAX_HEIGHT_FRC and CRELAX_HEIGHT_TYPE
+!     - LTRANS
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (module MODD_FRC)
+!!      
+!!
+!!    AUTHOR
+!!    ------
+!!	    Marc Georgelin Labo d'aerologie
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original 29/07/96 
+!!      29/07/96 (Pinty&Suhre) revised
+!!      18/11/96 J.-P. Pinty   addition of the translation
+!!      27/01/98 P. Bechtold   use tendency forcing
+!!                             add SST and surface pressure forcing
+!!      01/2004  V. Masson     surface externalization: removes SST forcing
+!!                   09/2017 Q.Rodier add LTEND_UV_FRC
+!!      03/2021 JL Redelsperger Parameters defining sfc forcing shape for idealized ocean case
+!!      06/2021 F. Couvreux    add LRELAX_UVMEAN_FRC
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+!USE MODD_TYPE_DATE
+!
+IMPLICIT NONE
+!
+!*            fields for FORCING
+!             ------------------
+!
+INTEGER,          SAVE                  :: NFRC     ! number of forcing profiles
+!TYPE (DATE_TIME), SAVE, DIMENSION(:), ALLOCATABLE :: TDTFRC ! date of
+                                                    !  each forcing profile
+!
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XUFRC,   &! geostrophic wind 
+					                       XVFRC,   &! components U and V
+					                       XWFRC     ! large scale vertical wind
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XTHFRC,  &! large scale TH profile
+					                       XRVFRC    ! large scale RV profile
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XGXTHFRC,&! large scale TH gradient
+                                           XGYTHFRC  ! along the X and Y axis
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XTENDTHFRC,&! large scale TH tendency
+                                           XTENDRVFRC  ! large scale RV tendency
+REAL, SAVE                              :: XUTRANS, &! horizontal components of
+                                           XVTRANS   !        a constant
+                                                     !   Galilean TRANSlation
+REAL, SAVE, DIMENSION(:), ALLOCATABLE   :: XPGROUNDFRC! surf. pressure 
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XTENDUFRC   ! large scale U tendency
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XTENDVFRC   ! large scale V tendency
+!
+!*            control parameters for FORCING
+!             ------------------------------
+!
+LOGICAL, SAVE     :: LGEOST_UV_FRC      ! enables geostrophic wind term
+LOGICAL, SAVE     :: LGEOST_TH_FRC      ! enables thermal wind advection
+LOGICAL, SAVE     :: LTEND_THRV_FRC     ! enables tendency forcing
+LOGICAL, SAVE     :: LTEND_UV_FRC       ! enables tendency forcing of the wind
+LOGICAL, SAVE     :: LVERT_MOTION_FRC   ! enables prescribed a forced vertical
+					                    ! transport for all prognostic variables
+LOGICAL, SAVE     :: LRELAX_THRV_FRC    ! enables temp. and humidity relaxation
+LOGICAL, SAVE     :: LRELAX_UV_FRC      ! enables  horizontal wind relaxation applied to the full wind field
+LOGICAL, SAVE     :: LRELAX_UVMEAN_FRC  ! enables  horizontal wind relaxation applied to the horiz. avg. wind
+!
+REAL,    SAVE     :: XRELAX_TIME_FRC    ! e-folding time for relaxation 
+REAL,    SAVE     :: XRELAX_HEIGHT_FRC  ! height below which relaxation
+                                        ! is never applied
+CHARACTER(len=4), SAVE :: CRELAX_HEIGHT_TYPE ! "THGR" relax. above maximal dTH/dz
+					                    ! (but always above XRELAX_HEIGHT_FRC)
+					                    ! "FIXE" relax. above XRELAX_HEIGHT_FRC
+!
+LOGICAL, SAVE     :: LTRANS             ! enables a Galilean translation of the
+                                        !         domain of simulation
+LOGICAL, SAVE     :: LPGROUND_FRC       ! enables surf. pressure forcing
+!
+LOGICAL, SAVE     :: LDEEPOC            ! activates sfc forcing for ideal ocean deep conv 
+REAL,    SAVE     :: XCENTX_OC          ! center of sfc forc for ideal ocean
+REAL,    SAVE     :: XRADX_OC           ! radius of sfc forc for ideal ocean
+REAL,    SAVE     :: XCENTY_OC          ! center of sfc forc for ideal ocean
+REAL,    SAVE     :: XRADY_OC           ! radius of sfc forc for ideal ocean
+!
+END MODULE MODD_FRC
diff --git a/src/arome/aux/modd_oceanh.f90 b/src/arome/aux/modd_oceanh.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6638752658492bb2a477b6781ca92fe885bfcf04
--- /dev/null
+++ b/src/arome/aux/modd_oceanh.f90
@@ -0,0 +1,47 @@
+!MNH_LIC Copyright 2021-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 MODD_OCEANH
+!     #################
+!
+!!****  *MODD_OCEAN* - declaration of variables used in ocean version
+!!
+!!    PURPOSE
+!!    -------
+!       Declarative module for the variables
+!!      at interface for OCEAN LES MESONH version including auto-coupling O-A LES
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    AUTHOR
+!!    ------
+!!    JL Redelsperger LOPS
+!!   
+!!    MODIFICATIONS
+!!    -------------
+!!      Original 03/2021
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+!USE MODD_TYPE_DATE
+!
+IMPLICIT NONE
+!
+!*            fields for Sea Sfc FORCINGs
+!             ------------------
+!
+INTEGER,          SAVE                  :: NFRCLT     ! number of sea surface forcings PLUS 1
+INTEGER,          SAVE                  :: NINFRT     ! Interval in second between forcings
+!TYPE (DATE_TIME), SAVE, DIMENSION(:), ALLOCATABLE :: TFRCLT ! date/time of sea surface forcings
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE ::  XSSUFL,XSSVFL,XSSTFL,XSSOLA ! Time evol Flux U V T Solar_Rad at sea surface
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XSSUFL_XY,XSSVFL_XY,XSSTFL_XY! XY flux shape
+REAL, SAVE, DIMENSION(:), ALLOCATABLE :: XSSUFL_T,XSSVFL_T,XSSTFL_T,XSSOLA_T ! given time forcing fluxes
+!
+END MODULE MODD_OCEANH
diff --git a/src/arome/aux/mode_gather_ll.F90 b/src/arome/aux/mode_gather_ll.F90
new file mode 100644
index 0000000000000000000000000000000000000000..dd922f659edf6f298379fa0a7acb29e7e3060acb
--- /dev/null
+++ b/src/arome/aux/mode_gather_ll.F90
@@ -0,0 +1,12 @@
+MODULE MODE_GATHER_ll
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE GATHERALL_FIELD_ll(HDIR,PSEND,PRECV,KRESP)
+CHARACTER(LEN=*),      INTENT(IN) :: HDIR
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PSEND
+REAL,DIMENSION(:,:,:), INTENT(INOUT):: PRECV
+INTEGER,               INTENT(INOUT):: KRESP
+
+CALL ABORT
+END SUBROUTINE GATHERALL_FIELD_ll
+END MODULE MODE_GATHER_ll
diff --git a/src/arome/aux/modi_second_mnh.F90 b/src/arome/aux/modi_second_mnh.F90
new file mode 100644
index 0000000000000000000000000000000000000000..aaa03bd45a5d4b51d58c4824b59e51184b62984c
--- /dev/null
+++ b/src/arome/aux/modi_second_mnh.F90
@@ -0,0 +1,7 @@
+MODULE MODI_SECOND_MNH
+INTERFACE
+SUBROUTINE SECOND_MNH(XT)
+REAL           :: XT
+END SUBROUTINE SECOND_MNH
+END INTERFACE
+END MODULE MODI_SECOND_MNH
diff --git a/src/arome/gmkpack_ignored_files b/src/arome/gmkpack_ignored_files
index 10c5e38ab406d69febb7865b66c5c4114f6bc608..1d3c455969a374ca79d8d836b2e5323c5fb1d728 100644
--- a/src/arome/gmkpack_ignored_files
+++ b/src/arome/gmkpack_ignored_files
@@ -85,3 +85,11 @@ phyex/micro/ini_budget.F90
 phyex/micro/modd_budget.F90
 phyex/micro/modi_budget.F90
 phyex/micro/modi_ini_budget.F90
+phyex/turb/tke_eps_sources.F90
+phyex/turb/turb_ver.F90
+phyex/turb/prandtl.F90
+phyex/turb/turb_ver_thermo_flux.F90
+phyex/turb/turb_ver_thermo_corr.F90
+phyex/turb/turb_ver_dyn_flux.F90
+phyex/turb/turb_ver_sv_flux.F90
+phyex/turb/turb_ver_sv_corr.F90
diff --git a/src/arome/turb/modd_turbn.F90 b/src/arome/turb/modd_turbn.F90
deleted file mode 100644
index 3d11a7b591941538033c61d4e9fa3e1f71b1d8d1..0000000000000000000000000000000000000000
--- a/src/arome/turb/modd_turbn.F90
+++ /dev/null
@@ -1,3 +0,0 @@
-MODULE MODD_TURB_n
- CHARACTER (LEN=4), SAVE  :: CTURBLEN='BL89'
-ENDMODULE MODD_TURB_n
diff --git a/src/arome/turb/modi_tke_eps_sources.F90 b/src/arome/turb/modi_tke_eps_sources.F90
deleted file mode 100644
index 28176fa63c0799e5609c458dc8fbeee2d669ebae..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_tke_eps_sources.F90
+++ /dev/null
@@ -1,60 +0,0 @@
-!     ######spl
-      MODULE MODI_TKE_EPS_SOURCES
-!     ###########################
-INTERFACE
-!
-      SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,  &
-                    & PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,        &
-                    & PTSTEP,PIMPL,PEXPL,                              &
-                    & HTURBLEN,HTURBDIM,                               &
-                    & TPFILE,OTURB_DIAG,                               &
-                    & PTP,PRTKES,PRTHLS,PCOEF_DISS,PTDIFF,PTDISS,      &
-                    & TBUDGETS, KBUDGETS,                              &
-                    & PEDR, PTR,PDISS,PRTKESM                          )
-                    !
-USE MODD_IO, ONLY: TFILEDATA
-USE MODD_BUDGET, ONLY : TBUDGETDATA
-!
-INTEGER,                 INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                 INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                 INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
-
-INTEGER,                 INTENT(IN)   ::  KMI          ! model index number  
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTKEM        ! TKE at t-deltat
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLM          ! mixing length         
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PLEPS        ! dissipative length
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRHODJ       ! density * grid volume
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                       ! metric coefficients
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PZZ          ! physical height w-pt
-REAL,                    INTENT(IN)   ::  PTSTEP       ! Time step 
-REAL,                    INTENT(IN)   ::  PEXPL, PIMPL ! Coef. temporal. disc.
-CHARACTER(LEN=4),        INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                       ! turbulence scheme
-CHARACTER(LEN=4),        INTENT(IN)   ::  HTURBLEN     ! kind of mixing length
-TYPE(TFILEDATA),         INTENT(IN)   ::  TPFILE       ! Output file
-LOGICAL,                 INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
-                                  ! diagnostic fields in the syncronous FM-file
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP          ! Dyn. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTRH
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTP          ! Ther. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTKES       ! RHOD * Jacobian *
-                                                       ! TKE at t+deltat
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTHLS       ! Source of Theta_l
-REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PCOEF_DISS   ! 1/(Cph*Exner)
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PTDIFF       ! Diffusion TKE term
-REAL, DIMENSION(:,:,:),  INTENT(OUT)  ::  PTDISS       ! Dissipation TKE term
-TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
-INTEGER, INTENT(IN) :: KBUDGETS
-REAL, DIMENSION(:,:,:),  INTENT(OUT), OPTIONAL  ::  PTR          ! Transport prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(OUT), OPTIONAL  ::  PDISS        ! Dissipation of TKE
-REAL, DIMENSION(:,:,:),  INTENT(OUT), OPTIONAL  ::  PEDR         ! EDR 
-REAL, DIMENSION(:,:,:),  INTENT(IN),  OPTIONAL  ::  PRTKESM      ! Advection source 
-!
-!
-!
-END SUBROUTINE TKE_EPS_SOURCES
-!
-END INTERFACE
-!
-END MODULE MODI_TKE_EPS_SOURCES
diff --git a/src/arome/turb/modi_turb_ver.F90 b/src/arome/turb/modi_turb_ver.F90
deleted file mode 100644
index 56dbe815d09ec06450d646e51345a6c120c157e6..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver.F90
+++ /dev/null
@@ -1,125 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER 
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER(KKA,KKU,KKL,KRR,KRRL,KRRI,                &
-                      OCLOSE_OUT,OTURB_FLX,                         &
-                      HTURBDIM,HTOM,PIMPL,PEXPL,                    & 
-                      PTSTEP_UVW,PTSTEP_MET, PTSTEP_SV,             &
-                      HFMFILE,HLUOUT,                               &
-                      PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
-                      PCOSSLOPE,PSINSLOPE,                          &
-                      PRHODJ,PTHVREF,                               &
-                      PSFTHM,PSFRM,PSFSVM,PSFTHP,PSFRP,PSFSVP,      &
-                      PCDUEFF,PTAU11M,PTAU12M,PTAU33M,              &
-                      PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
-                      PTKEM,PLM,PLENGTHM,PLENGTHH,PLEPS,MFMOIST,    &
-                      PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE,    &
-                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,PBL_DEPTH,        &
-                      PSBL_DEPTH,PLMO,                              &
-                      PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,             &
-                      PDP,PTP,PSIGS,PWTH,PWRC,PWSV          )
-
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=ARO
-INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
-INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
-INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the 
-                                                      ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-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, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX, PDYY, PDZZ, PDZX, PDZY 
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PZZ          ! altitudes at flux points
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCOSSLOPE    ! cosinus of the angle 
-                                      ! between i and the slope vector
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE    ! sinus of the angle 
-                                      ! between i and the slope vector
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
-
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
-                                                      ! Potential Temperature 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHM,PSFRM ! surface fluxes at time
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSFSVM       ! t - deltat 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHP,PSFRP ! surface fluxes at time
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSFSVP       ! t + deltat 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCDUEFF     ! Cd * || u || at time t
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU11M      ! <uu> in the axes linked 
-       ! to the maximum slope direction and the surface normal and the binormal 
-       ! at time t - dt
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU12M      ! <uv> in the same axes
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU33M      ! <ww> in the same axes
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PUM,PVM,PWM,PTHLM 
-  ! Wind and potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios 
-                                                      ! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! scalar var. at t-Delta t
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PUSLOPEM     ! wind component along the 
-                                     ! maximum slope direction
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PVSLOPEM     ! wind component along the 
-                                     ! direction normal to the maximum slope one
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLENGTHM     ! Turb. mixing length momentum   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLENGTHH     ! Turb. mixing length heat/moisture  
-
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM        ! normalized 
-                  ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFRAC_ICE    ! ri fraction of rc+ri
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWTH        ! d(w'2th' )/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWR         ! d(w'2r'  )/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTH2        ! d(w'th'2 )/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFR2         ! d(w'r'2  )/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTHR        ! d(w'th'r')/dz
-REAL, DIMENSION(:,:),   INTENT(INOUT)::  PBL_DEPTH    ! BL depth
-REAL, DIMENSION(:,:),   INTENT(INOUT)::  PSBL_DEPTH   ! SBL depth
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PLMO         ! Monin-Obukhov length
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT)   ::  PRUS, PRVS, PRWS, PRTHLS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRSVS,PRRS
-                            ! cumulated sources for the prognostic variables
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PDP,PTP   ! Dynamic and thermal
-                                                   ! TKE production terms
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PSIGS     ! Vert. part of Sigma_s at t
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTH       ! heat flux
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWRC       ! cloud water flux
-REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
-
-!
-!
-END SUBROUTINE TURB_VER
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER
diff --git a/src/arome/turb/modi_turb_ver_dyn_flux.F90 b/src/arome/turb/modi_turb_ver_dyn_flux.F90
deleted file mode 100644
index a58d4f8a34cbd0d308c9b370765f370458e33e1c..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver_dyn_flux.F90
+++ /dev/null
@@ -1,86 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER_DYN_FLUX 
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL,                     &
-                      OCLOSE_OUT,OTURB_FLX,KRR,                     &
-                      HTURBDIM,PIMPL,PEXPL,                         &
-                      PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
-                      PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
-                      PCOSSLOPE,PSINSLOPE,                          &
-                      PRHODJ,                                       &
-                      PCDUEFF,PTAU11M,PTAU12M,PTAU33M,              &
-                      PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
-                      PTKEM,PLM,MFMOIST,PWU,PWV,                    &
-                      PRUS,PRVS,PRWS,                               &
-                      PDP,PTP                                )
-!
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=AR 
-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
-INTEGER,                INTENT(IN)   ::  KRR          ! number of moist var.
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the 
-                                                      ! turbulence scheme
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-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)   ::  PZZ          ! altitude of flux points
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCOSSLOPE    ! cosinus of the angle 
-                                      ! between i and the slope vector
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE    ! sinus of the angle 
-                                      ! between i and the slope vector
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
-
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCDUEFF     ! Cd * || u || at time t
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU11M      ! <uu> in the axes linked 
-       ! to the maximum slope direction and the surface normal and the binormal 
-       ! at time t - dt
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU12M      ! <uv> in the same axes
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU33M      ! <ww> in the same axes
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PUM,PVM,PWM,PTHLM
-  ! Wind at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PUSLOPEM     ! wind component along the 
-                                     ! maximum slope direction
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PVSLOPEM     ! wind component along the 
-                                     ! direction normal to the maximum slope one
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PWU          ! momentum flux u'w'
-REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PWV          ! momentum flux v'w'
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT)   ::  PRUS, PRVS, PRWS
-                            ! cumulated sources for the prognostic variables
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT)::  PDP,PTP   ! Dynamic and thermal
-                                                   ! TKE production terms
-!
-!
-!
-END SUBROUTINE TURB_VER_DYN_FLUX
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_DYN_FLUX
diff --git a/src/arome/turb/modi_turb_ver_sv_corr.F90 b/src/arome/turb/modi_turb_ver_sv_corr.F90
deleted file mode 100644
index 31a993729572ee1f9f42393947976daf65d5a093..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver_sv_corr.F90
+++ /dev/null
@@ -1,44 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER_SV_CORR
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,                    &
-                      PDZZ,                                         &
-                      PTHLM,PRM,PTHVREF,                            &
-                      PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PPHI3,PPSI3,  &
-                      PWM,PSVM,                                     &
-                      PTKEM,PLM,PLEPS,PPSI_SV                       )
-!
-INTEGER,                 INTENT(IN)  ::  KKA, KKU ! near ground and uppest atmosphere array indexes
-INTEGER,                 INTENT(IN)  ::  KKL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-INTEGER,                INTENT(IN)   ::  KRR          ! number of moist var.
-INTEGER,                INTENT(IN)   ::  KRRL         ! number of liquid var.
-INTEGER,                INTENT(IN)   ::  KRRI         ! number of ice var.
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDZZ 
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM        ! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios at t-Delta t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! reference Thv
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM        ! normalized 
-                  ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPHI3        ! Inv.Turb.Sch.for temperature
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPSI3        ! Inv.Turb.Sch.for humidity
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PWM          ! w at time t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! scalar var. at t-Delta t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PPSI_SV      ! Inv.Turb.Sch.for scalars
-                            ! cumulated sources for the prognostic variables
-!
-!
-END SUBROUTINE TURB_VER_SV_CORR
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_SV_CORR
diff --git a/src/arome/turb/modi_turb_ver_sv_flux.F90 b/src/arome/turb/modi_turb_ver_sv_flux.F90
deleted file mode 100644
index dd10b563a2ff8686a345e31a6c75f662a56a86b2..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver_sv_flux.F90
+++ /dev/null
@@ -1,65 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER_SV_FLUX 
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL,                      &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,                &
-                      PIMPL,PEXPL,                                  &
-                      PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
-                      PDZZ,PDIRCOSZW,                               &
-                      PRHODJ,PWM,                                   &
-                      PSFSVM,PSFSVP,                                &
-                      PSVM,                                         &
-                      PTKEM,PLM,MFMOIST,PPSI_SV,                    &
-                      PRSVS,PWSV                                    )
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=AR
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                      ! turbulence scheme
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-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)   ::  PDZZ 
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
-
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSFSVM       ! t - deltat 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSFSVP       ! t + deltat 
-!
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! scalar var. at t-Delta t
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PWM          ! vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PPSI_SV      ! Inv.Turb.Sch.for scalars
-!
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) ::  PRSVS
-                            ! cumulated sources for the prognostic variables
-REAL, DIMENSION(:,:,:,:), INTENT(OUT)  :: PWSV        ! scalar flux
-
-!
-!
-END SUBROUTINE TURB_VER_SV_FLUX
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_SV_FLUX
diff --git a/src/arome/turb/modi_turb_ver_thermo_corr.F90 b/src/arome/turb/modi_turb_ver_thermo_corr.F90
deleted file mode 100644
index 43d35b7707fc90cc4dee3546333e1cbaff4ef77c..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver_thermo_corr.F90
+++ /dev/null
@@ -1,107 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER_THERMO_CORR
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
-                      PIMPL,PEXPL,                                  &
-                      HFMFILE,HLUOUT,                               &
-                      PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,           &
-                      PRHODJ,PTHVREF,                               &
-                      PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
-                      PWM,PTHLM,PRM,PSVM,                           &
-                      PTKEM,PLM,PLEPS,                              &
-                      PLOCPEXNM,PATHETA,PAMOIST,PSRCM,              &
-                      PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3,  &
-                      PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
-                      PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
-                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
-                      PTHLP,PRP,MFMOIST,PSIGS                  )
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=AR 
-INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
-INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
-INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                      ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-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)   ::  PDZZ, PDXX, PDYY, PDZX, PDZY
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
-
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
-                                                      ! Potential Temperature 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHM,PSFRM ! surface fluxes at time
-!                                                     ! t - deltat 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHP,PSFRP ! surface fluxes at time
-!                                                     ! t + deltat 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PWM 
-! Vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM 
-! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios 
-                                                      ! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! Mixing ratios 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM        ! normalized 
-                  ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBETA        ! buoyancy coefficient
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSQRT_TKE    ! sqrt(e)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDTH_DZ      ! d(th)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDR_DZ       ! d(rt)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2TH3     ! 3D Redeslperger number R*2_th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2R3      ! 3D Redeslperger number R*2_r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2THR3    ! 3D Redeslperger number R*2_thr
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBLL_O_E     ! beta * Lk * Leps / tke
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PETHETA      ! Coefficient for theta in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PEMOIST      ! Coefficient for r in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDTH1      ! 1D Redelsperger number for Th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDR1       ! 1D Redelsperger number for r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPHI3        ! Prandtl number for temperature
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPSI3        ! Prandtl number for vapor
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PD           ! Denominator in Prandtl numbers
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWTH        ! d(w'2th' )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWR         ! d(w'2r'  )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTH2        ! d(w'th'2 )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFR2         ! d(w'r'2  )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTHR        ! d(w'th'r')/dz (at mass point)
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHLP      ! guess of thl at t+ deltat
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRP        ! guess of r at t+ deltat
-!
-REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PSIGS     ! Vert. part of Sigma_s at t
-!
-!
-!
-END SUBROUTINE TURB_VER_THERMO_CORR
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_THERMO_CORR
diff --git a/src/arome/turb/modi_turb_ver_thermo_flux.F90 b/src/arome/turb/modi_turb_ver_thermo_flux.F90
deleted file mode 100644
index 19a050e790bf82229cd2acd0a5b3a670045e83c5..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_turb_ver_thermo_flux.F90
+++ /dev/null
@@ -1,118 +0,0 @@
-!     ######spl
-     MODULE MODI_TURB_VER_THERMO_FLUX 
-!    ####################
-!
-INTERFACE 
-!
-      SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
-                      PIMPL,PEXPL,                                  &
-                      PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
-                      PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
-                      PRHODJ,PTHVREF,                               &
-                      PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
-                      PWM,PTHLM,PRM,PSVM,                           &
-                      PTKEM,PLM,PLEPS,                              &
-                      PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE,    &
-                      PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3,  &
-                      PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
-                      PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
-                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,MFMOIST,PBL_DEPTH,&
-                      PWTHV,PRTHLS,PRRS,PTHLP,PRP,PTP,PWTH,PWRC )
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=AR O
-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.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                      ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-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)   ::  PDZZ, PDXX, PDYY, PDZX, PDZY
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PZZ          ! altitudes
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
-
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
-                                                      ! Potential Temperature 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHM,PSFRM ! surface fluxes at time
-!                                                     ! t - deltat 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHP,PSFRP ! surface fluxes at time
-!                                                     ! t + deltat 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PWM 
-! Vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM 
-! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios 
-                                                      ! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! Mixing ratios 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM        ! normalized 
-                  ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFRAC_ICE    ! ri fraction of rc+ri
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBETA        ! buoyancy coefficient
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSQRT_TKE    ! sqrt(e)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDTH_DZ      ! d(th)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDR_DZ       ! d(rt)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2TH3     ! 3D Redeslperger number R*2_th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2R3      ! 3D Redeslperger number R*2_r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2THR3    ! 3D Redeslperger number R*2_thr
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBLL_O_E     ! beta * Lk * Leps / tke
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PETHETA      ! Coefficient for theta in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PEMOIST      ! Coefficient for r in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDTH1      ! 1D Redelsperger number for Th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDR1       ! 1D Redelsperger number for r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPHI3        ! Prandtl number for temperature
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPSI3        ! Prandtl number for vapor
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PD           ! Denominator in Prandtl numbers
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWTH        ! d(w'2th' )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWR         ! d(w'2r'  )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTH2        ! d(w'th'2 )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFR2         ! d(w'r'2  )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTHR        ! d(w'th'r')/dz (at mass point)
-REAL, DIMENSION(:,:),   INTENT(INOUT)::  PBL_DEPTH    ! BL depth
-REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTHV         ! buoyancy flux
-!
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHLS     ! cumulated source for theta
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS       ! cumulated source for rt
-REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTHLP      ! guess of thl at t+ deltat
-REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRP        ! guess of r at t+ deltat
-!
-REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PTP       ! Dynamic and thermal
-                                                     ! TKE production terms
-!
-REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PWTH       ! heat flux
-REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PWRC       ! cloud water flux
-!
-!
-END SUBROUTINE TURB_VER_THERMO_FLUX
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_THERMO_FLUX
diff --git a/src/common/aux/modd_blowsnow.f90 b/src/common/aux/modd_blowsnow.f90
new file mode 100644
index 0000000000000000000000000000000000000000..36180eadd60e261aa7b39b19e96edbda940c31f6
--- /dev/null
+++ b/src/common/aux/modd_blowsnow.f90
@@ -0,0 +1,80 @@
+!MNH_LIC Copyright 1994-2018 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 MODD_BLOWSNOW
+!!     ######################
+!!
+!!     PURPOSE
+!!     -------
+!!
+!!     Declaration of variables and types for the blowing snow scheme
+!!
+!!     METHOD
+!!     ------
+!!
+!!
+!!     REFERENCE
+!!     ---------
+!!     Etudes du transport de la neige par le vent en conditions alpines : 
+!!     Observations et simulations à l'aide d'un modèle couplé atmosphère/
+!!     manteau neigeux (Thèse, Uni. Paris Est, 2012)
+!!
+!!
+!!     AUTHOR
+!!     ------
+!!     Vincent Vionnet (CNRM)
+!!
+!!
+!!     MODIFICATIONS
+!!     -------------
+!!
+!!--------------------------------------------------------------------
+!!     DECLARATIONS
+!!     ------------
+IMPLICIT NONE
+
+LOGICAL      :: LBLOWSNOW  = .FALSE.   ! switch to active pronostic blowing snow
+!
+INTEGER      :: NBLOWSNOW3D = 2 ! Number of blowing snow variables
+! as scalar in Meso-NH. The curent version of the model use two scalars:
+!				- Number concentration (#/kg)
+!				- Mass concentration (kg/kg)
+
+INTEGER     :: NBLOWSNOW_2D = 3  ! Number of 2D blowing snow variables
+! adected in Meso-NH. The curent version of the model advectes three variables:
+!             - total number concentration in Canopy
+!             - total mass concentration in Canopy
+!             - equivalent concentration in the saltation layer
+!
+REAL            :: XALPHA_SNOW ! Gamma distribution shape factor
+!
+REAL            :: XRSNOW  ! Ratio between diffusion coefficient for scalar
+                           ! variables and blowing snow variables
+                           ! RSNOW = KSCA/KSNOW = 4. (if Redelsperger-Sommeria (1981) used in ini_cturb)
+                           ! RSNOW = KSCA/KSNOW = 2.5 ( if Cheng-Canuto-Howard (2002) used in ini_cturb)
+                           ! Cheng-Canuto-Howard (2002) is the default in MNH V5.3
+                           ! See Vionnet (PhD, 2012, In French) and Vionnet et al (TC, 2014)
+                           ! for a complete dicsussion
+CHARACTER(LEN=6),DIMENSION(:),ALLOCATABLE  :: CSNOWNAMES
+
+CHARACTER(LEN=6),DIMENSION(2), PARAMETER  :: YPSNOW_INI = &
+     (/'SNWM01','SNWM02'/)
+!
+CHARACTER(LEN=6),DIMENSION(3), PARAMETER  :: YPBLOWSNOW_2D = &
+     (/'SNWCNU','SNWCMA','SNWCSA' /)
+
+CHARACTER(LEN=4)  :: CSNOWSEDIM ! type of formulation for snow
+!              sedimentation : MITC : Mitchell (1996)
+!                              CARR : Carrier's drag coefficient (cf PIEKTUK)
+!                              TABC : Tabulated values from Carrier's drag coefficient
+!                              NONE : no seidmentation
+!Minimal mean radius (um) 
+REAL           :: XINIRADIUS_SNW = 5.e-6
+!Minimum allowed number concentration (#/m3)
+REAL           :: XN0MIN_SNW    =  1
+!
+!
+END MODULE MODD_BLOWSNOW
diff --git a/src/common/aux/modd_dimn.f90 b/src/common/aux/modd_dimn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..622541f68f37ad130a45ecd4169932eae50f2967
--- /dev/null
+++ b/src/common/aux/modd_dimn.f90
@@ -0,0 +1,87 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 modd 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+!     ##################
+      MODULE MODD_DIM_n
+!     ##################
+!
+!!****  *MODD_DIM$n* - declaration of dimensions
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to specify  the dimensions 
+!     of the data arrays.   
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (module MODD_DIMn)
+!!      Technical Specifications Report of the Meso-NH (chapters 2 and 3)
+!!          
+!!    AUTHOR
+!!    ------
+!!	V. Ducrocq   *Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/05/94     
+!!      Modifications 13/08/98 (V. Ducrocq) // NIINF .. NJSUP are no more used in the init part                
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+IMPLICIT NONE
+
+TYPE DIM_t
+  INTEGER :: NIMAX,NJMAX,NKMAX  !  Dimensions respectively  in x , 
+                              ! y ,  z directions of the physical sub-domain.
+  INTEGER :: NIMAX_ll,NJMAX_ll  !  Dimensions respectively  in x and y
+                                   ! directions of the physical domain
+  INTEGER :: NIINF, NISUP       !  Lower bound and upper bound of the arrays 
+                                   ! in x direction 
+  INTEGER :: NJINF, NJSUP       !  Lower bound and upper bound of the arrays 
+                                   ! in y direction
+!
+END TYPE DIM_t
+
+TYPE(DIM_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: DIM_MODEL
+
+INTEGER, POINTER :: NIMAX=>NULL(),NJMAX=>NULL(),NKMAX=>NULL()
+INTEGER, POINTER :: NIMAX_ll=>NULL(),NJMAX_ll=>NULL()
+INTEGER, POINTER :: NIINF=>NULL(), NISUP=>NULL()
+INTEGER, POINTER :: NJINF=>NULL(), NJSUP=>NULL()
+
+CONTAINS
+
+SUBROUTINE DIM_GOTO_MODEL(KFROM, KTO)
+INTEGER, INTENT(IN) :: KFROM, KTO
+!
+! Save current state for allocated arrays
+!
+! Current model is set to model KTO
+NIMAX=>DIM_MODEL(KTO)%NIMAX
+NJMAX=>DIM_MODEL(KTO)%NJMAX
+NKMAX=>DIM_MODEL(KTO)%NKMAX
+NIMAX_ll=>DIM_MODEL(KTO)%NIMAX_ll
+NJMAX_ll=>DIM_MODEL(KTO)%NJMAX_ll
+NIINF=>DIM_MODEL(KTO)%NIINF
+NISUP=>DIM_MODEL(KTO)%NISUP
+NJINF=>DIM_MODEL(KTO)%NJINF
+NJSUP=>DIM_MODEL(KTO)%NJSUP
+
+END SUBROUTINE DIM_GOTO_MODEL
+
+END MODULE MODD_DIM_n
diff --git a/src/common/aux/modd_gridn.f90 b/src/common/aux/modd_gridn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..055d3c88f76b4634b987607db83365a12e58710f
--- /dev/null
+++ b/src/common/aux/modd_gridn.f90
@@ -0,0 +1,67 @@
+!MNH_LIC Copyright 1994-2018 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 MODD_GRID_n
+!     ##################
+!
+!!****  *MODD_GRID$n* - declaration of grid variables
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to declare  the variables
+!     describing the grid. 
+!    
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (module MODD_GRIDn)
+!!      Technical Specifications Report of the Meso-NH (chapters 2 and 3)
+!!
+!!    AUTHOR
+!!    ------
+!!	V. Ducrocq   *Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/05/94                      
+!!      J. Stein    15/11/95  add the slope angle
+!!      V. Ducrocq   13/08/98  // : add XLATOR_ll and XLONOR_ll       
+!!      V. Masson   nov 2004  supress XLATOR,XLONOR,XLATOR_ll,XLONOR_ll
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+IMPLICIT NONE
+
+REAL, DIMENSION(:,:),  POINTER :: XLON=>NULL(),XLAT=>NULL() ! Longitude and latitude  
+REAL, DIMENSION(:),    POINTER :: XXHAT=>NULL()             ! Position x in the conformal or cartesian plane
+REAL, DIMENSION(:),    POINTER :: XYHAT=>NULL()             ! Position y in the conformal or cartesian plane
+REAL, DIMENSION(:),    POINTER :: XDXHAT=>NULL()            ! horizontal stretching in x
+REAL, DIMENSION(:),    POINTER :: XDYHAT=>NULL()            ! horizontal stretching in y
+REAL, DIMENSION(:,:),  POINTER :: XMAP=>NULL()              ! Map factor 
+REAL, DIMENSION(:,:),  POINTER :: XZS=>NULL()               ! orography
+REAL, DIMENSION(:,:,:),POINTER :: XZZ=>NULL()               ! height z 
+REAL,                  POINTER :: XZTOP=>NULL()             ! model top (m)
+REAL, DIMENSION(:),    POINTER :: XZHAT=>NULL()             ! height level without orography
+REAL, DIMENSION(:,:),  POINTER :: XDIRCOSXW=>NULL(),XDIRCOSYW=>NULL(),XDIRCOSZW=>NULL() ! director cosinus of the normal
+                                                                                        ! to the ground surface
+REAL, DIMENSION(:,:),  POINTER  :: XCOSSLOPE=>NULL()         ! cosinus of the angle between i and the slope vector
+REAL, DIMENSION(:,:),  POINTER  :: XSINSLOPE=>NULL()         ! sinus   of the angle between i and the slope vector
+! Quantities for SLEVE vertical coordinate
+LOGICAL,               POINTER  :: LSLEVE=>NULL()            ! Logical for SLEVE coordinate
+REAL,                  POINTER  :: XLEN1=>NULL()             ! Decay scale for smooth topography
+REAL,                  POINTER  :: XLEN2=>NULL()             ! Decay scale for small-scale topography deviation
+REAL, DIMENSION(:,:),  POINTER  :: XZSMT=>NULL()             ! smooth orography for SLEVE coordinate
+
+END MODULE MODD_GRID_n
diff --git a/src/common/aux/modd_metricsn.f90 b/src/common/aux/modd_metricsn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..33cec104e0fcc1f51650c617136afecbb7f8edd0
--- /dev/null
+++ b/src/common/aux/modd_metricsn.f90
@@ -0,0 +1,78 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 modd 2006/06/27 14:20:29
+!-----------------------------------------------------------------
+!     #####################
+      MODULE MODD_METRICS_n
+!     #####################
+!
+!!****  *MODD_METRICS$n* - metric coefficients
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to declare the 
+!     metric coefficients. 
+!    
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!
+!!    AUTHOR
+!!    ------
+!!	P. Jabouille   *Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/04/99                      
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+IMPLICIT NONE
+
+TYPE METRICS_t
+  REAL, DIMENSION(:,:,:), POINTER :: XDXX=>NULL(),XDZX=>NULL(), &
+                                  XDYY=>NULL(),XDZY=>NULL(),XDZZ=>NULL()
+                                            !metric coefficients
+END TYPE METRICS_t
+
+TYPE(METRICS_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: METRICS_MODEL
+
+REAL, DIMENSION(:,:,:), POINTER :: XDXX=>NULL(),XDZX=>NULL(), &
+                                  XDYY=>NULL(),XDZY=>NULL(),XDZZ=>NULL()
+
+CONTAINS
+
+SUBROUTINE METRICS_GOTO_MODEL(KFROM, KTO)
+INTEGER, INTENT(IN) :: KFROM, KTO
+!
+! Save current state for allocated arrays
+METRICS_MODEL(KFROM)%XDXX=>XDXX
+METRICS_MODEL(KFROM)%XDZX=>XDZX
+METRICS_MODEL(KFROM)%XDYY=>XDYY
+METRICS_MODEL(KFROM)%XDZY=>XDZY
+METRICS_MODEL(KFROM)%XDZZ=>XDZZ
+!
+! Current model is set to model KTO
+XDXX=>METRICS_MODEL(KTO)%XDXX
+XDZX=>METRICS_MODEL(KTO)%XDZX
+XDYY=>METRICS_MODEL(KTO)%XDYY
+XDZY=>METRICS_MODEL(KTO)%XDZY
+XDZZ=>METRICS_MODEL(KTO)%XDZZ
+
+END SUBROUTINE METRICS_GOTO_MODEL
+
+END MODULE MODD_METRICS_n
diff --git a/src/common/aux/modd_ref.f90 b/src/common/aux/modd_ref.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3d1e6025ba5ac880c2d388e7adfe3be615ad2099
--- /dev/null
+++ b/src/common/aux/modd_ref.f90
@@ -0,0 +1,58 @@
+!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 MODD_REF
+!     ###############
+!
+!!****  *MODD_REF* - declaration of reference state  profile
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to declare  the vertical
+!     profile of  the reference state, used for the anelastic 
+!     approximation. 
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (module MODD_REF)
+!!      Technical Specifications Report of the Meso-NH (chapters 2 and 3)
+!!
+!!    AUTHOR
+!!    ------
+!!	V. Ducrocq   *Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original   07/06/94   
+!!                    07/13 (C.Lac) Add LBOUSS
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+IMPLICIT NONE
+!          
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XRHODREFZ ! rhod(z) for reference
+                                             ! state without orography
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XTHVREFZ  ! Thetav(z) for reference
+                                             ! state without orography    
+REAL,SAVE                            :: XEXNTOP   ! Exner function at model top 
+!
+! For coupled A-O case
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XRHODREFZO! rhod(z) for ocean ref state in coupled mode
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XTHVREFZO !Thetav(z) for ocean ref state in coupled mode
+REAL,SAVE                            :: XEXNTOPO   ! Exner function at ocean  model top in coupled mode
+!
+LOGICAL, SAVE                        :: LBOUSS    ! Boussinesq approximation
+LOGICAL, SAVE   ::LCOUPLES ! AUTOCOUPLED ATMS-OCEAN LES VERSION
+! 
+END MODULE MODD_REF
diff --git a/src/common/aux/modd_turbn.f90 b/src/common/aux/modd_turbn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..8c35fd9d4be7bd61167ee2e6aba4a4fb40e521a5
--- /dev/null
+++ b/src/common/aux/modd_turbn.f90
@@ -0,0 +1,211 @@
+!MNH_LIC Copyright 1995-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 MODD_TURB_n
+!     ##################
+!
+!!****  *MODD_TURB$n* - declaration of turbulence scheme free parameters
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to declare the
+!     variables that may be set by namelist for the turbulence scheme
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (module MODD_PARAMn)
+!!          
+!!    AUTHOR
+!!    ------
+!!	    J. Cuxart and J. Stein       * I.N.M. and Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    January 9, 1995                   
+!!      J.Cuxart    February 15, 1995 add the switches for diagnostic storages
+!!      J.M. Carriere May  15, 1995 add the subgrid condensation
+!!      M. Tomasini Jul  05, 2001 add the subgrid autoconversion
+!!      P. Bechtold Feb 11, 2002    add switch for Sigma_s computation
+!!      P. Jabouille Apr 4, 2002    add switch for Sigma_s convection
+!!      V. Masson    Nov 13 2002    add switch for SBL lengths
+!!                   May   2006    Remove KEPS
+!!      C.Lac        Nov 2014      add terms of TKE production for LES diag
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!      D. Ricard     May 2021      add the switches for Leonard terms
+!!    JL Redelsperger  03/2021   Add O-A flux for auto-coupled LES case
+!!
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+IMPLICIT NONE
+
+TYPE TURB_t
+! 
+! 
+  REAL               :: XIMPL     ! implicitness degree for the vertical terms of 
+                                     ! the turbulence scheme
+  REAL               :: XKEMIN      ! mimimum value for the TKE                                  
+  REAL               :: XCEDIS      ! Constant for dissipation of Tke                            
+  REAL               :: XCADAP      ! Coefficient for ADAPtative mixing length
+  CHARACTER (LEN=4)  :: CTURBLEN  ! type of length used for the closure
+                                     ! 'BL89' Bougeault and Lacarrere scheme
+                                     ! 'DELT' length = ( volum) ** 1/3
+  CHARACTER (LEN=4)  :: CTURBDIM  ! dimensionality of the turbulence scheme
+                                     ! '1DIM' for purely vertical computations
+                                     ! '3DIM' for computations in the 3 
+                                     ! directions
+  LOGICAL            :: LTURB_FLX ! logical switch for the storage of all  
+                                     ! the turbulent fluxes
+  LOGICAL            :: LTURB_DIAG! logical switch for the storage of some 
+                                     ! turbulence related diagnostics
+  LOGICAL            :: LSUBG_COND! Switch for subgrid condensation 
+  LOGICAL            :: LSIGMAS   ! Switch for using Sigma_s from turbulence scheme
+  LOGICAL            :: LSIG_CONV ! Switch for computing Sigma_s due to convection
+!
+  LOGICAL            :: LRMC01    ! Switch for computing separate mixing
+!                                    ! and dissipative length in the SBL
+!                                    ! according to Redelsperger, Mahe &
+!                                    ! Carlotti 2001
+  CHARACTER(LEN=4)   :: CTOM      ! type of Third Order Moments
+                                  ! 'NONE' none
+                                  ! 'TM06' Tomas Masson 2006
+  CHARACTER(LEN=4)   :: CSUBG_AUCV ! type of subgrid rc->rr autoconv. method
+  CHARACTER(LEN=80)  :: CSUBG_AUCV_RI ! type of subgrid ri->rs autoconv. method
+  CHARACTER(LEN=80)  :: CCONDENS ! subrgrid condensation PDF
+  CHARACTER(LEN=4)   :: CLAMBDA3 ! lambda3 choice for subgrid cloud scheme
+  CHARACTER(LEN=80)  :: CSUBG_MF_PDF ! PDF to use for MF cloud autoconversions
+
+!  REAL, DIMENSION(:,:), POINTER :: XBL_DEPTH=>NULL() ! BL depth for TOMS computations
+!  REAL, DIMENSION(:,:), POINTER :: XSBL_DEPTH=>NULL()! SurfaceBL depth for RMC01 computations
+!  REAL, DIMENSION(:,:,:), POINTER :: XWTHVMF=>NULL()! Mass Flux vert. transport of buoyancy
+  REAL               :: VSIGQSAT  ! coeff applied to qsat variance contribution
+  REAL, DIMENSION(:,:,:), POINTER :: XDYP=>NULL()    ! Dynamical production of Kinetic energy
+  REAL, DIMENSION(:,:,:), POINTER :: XTHP=>NULL()    ! Thermal production of Kinetic energy
+  REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL()    ! Transport production of Kinetic energy
+  REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL()    ! Dissipation of Kinetic energy
+  REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL()    ! Mixing length
+  REAL, DIMENSION(:,:,:), POINTER :: XSSUFL_C=>NULL() ! O-A interface flux for u
+  REAL, DIMENSION(:,:,:), POINTER :: XSSVFL_C=>NULL() ! O-A interface flux for v
+  REAL, DIMENSION(:,:,:), POINTER :: XSSTFL_C=>NULL() ! O-A interface flux for theta
+  REAL, DIMENSION(:,:,:), POINTER :: XSSRFL_C=>NULL() ! O-A interface flux for vapor
+  LOGICAL            :: LHGRAD ! logical switch for the computation of the Leornard Terms
+  REAL               :: XCOEFHGRADTHL  ! coeff applied to thl contribution
+  REAL               :: XCOEFHGRADRM  ! coeff applied to mixing ratio contribution
+  REAL               :: XALTHGRAD  ! altitude from which to apply the Leonard terms
+  REAL               :: XCLDTHOLD  ! cloud threshold to apply the Leonard terms
+                                   ! negative value : applied everywhere
+                                   ! 0.000001 applied only inside the clouds ri+rc > 10**-6 kg/kg
+!
+END TYPE TURB_t
+
+TYPE(TURB_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: TURB_MODEL
+
+REAL, POINTER :: XIMPL=>NULL()
+REAL, POINTER :: XKEMIN=>NULL()
+REAL, POINTER :: XCEDIS=>NULL()
+REAL, POINTER :: XCADAP=>NULL()
+CHARACTER (LEN=4), POINTER :: CTURBLEN=>NULL()
+CHARACTER (LEN=4), POINTER :: CTURBDIM=>NULL()
+LOGICAL, POINTER :: LTURB_FLX=>NULL()
+LOGICAL, POINTER :: LTURB_DIAG=>NULL()
+LOGICAL, POINTER :: LSUBG_COND=>NULL()
+LOGICAL, POINTER :: LSIGMAS=>NULL()
+LOGICAL, POINTER :: LSIG_CONV=>NULL()
+LOGICAL, POINTER :: LRMC01=>NULL()
+CHARACTER(LEN=4),POINTER :: CTOM=>NULL()
+CHARACTER(LEN=4),POINTER :: CSUBG_AUCV=>NULL()
+CHARACTER(LEN=80),POINTER :: CSUBG_AUCV_RI=>NULL()
+CHARACTER(LEN=80),POINTER :: CCONDENS=>NULL()
+CHARACTER(LEN=4),POINTER :: CLAMBDA3=>NULL()
+CHARACTER(LEN=80),POINTER :: CSUBG_MF_PDF=>NULL()
+REAL, DIMENSION(:,:), POINTER :: XBL_DEPTH=>NULL()
+REAL, DIMENSION(:,:), POINTER :: XSBL_DEPTH=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XWTHVMF=>NULL()
+REAL, POINTER :: VSIGQSAT=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XDYP=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XTHP=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSUFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSVFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSTFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSRFL_C=>NULL()
+LOGICAL, POINTER :: LHGRAD=>NULL()
+REAL, POINTER :: XCOEFHGRADTHL=>NULL()
+REAL, POINTER :: XCOEFHGRADRM=>NULL()
+REAL, POINTER :: XALTHGRAD=>NULL()
+REAL, POINTER :: XCLDTHOLD=>NULL()
+
+CONTAINS
+
+SUBROUTINE TURB_GOTO_MODEL(KFROM, KTO)
+INTEGER, INTENT(IN) :: KFROM, KTO
+!
+! Save current state for allocated arrays
+!
+!TURB_MODEL(KFROM)%XBL_DEPTH=>XBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
+!TURB_MODEL(KFROM)%XSBL_DEPTH=>XSBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
+!TURB_MODEL(KFROM)%XWTHVMF=>XWTHVMF !Done in FIELDLIST_GOTO_MODEL
+TURB_MODEL(KFROM)%XDYP=>XDYP 
+TURB_MODEL(KFROM)%XTHP=>XTHP 
+TURB_MODEL(KFROM)%XTR=>XTR 
+TURB_MODEL(KFROM)%XDISS=>XDISS
+TURB_MODEL(KFROM)%XLEM=>XLEM
+TURB_MODEL(KFROM)%XSSUFL_C=>XSSUFL_C
+TURB_MODEL(KFROM)%XSSVFL_C=>XSSVFL_C
+TURB_MODEL(KFROM)%XSSTFL_C=>XSSTFL_C
+TURB_MODEL(KFROM)%XSSRFL_C=>XSSRFL_C
+!
+! Current model is set to model KTO
+XIMPL=>TURB_MODEL(KTO)%XIMPL
+XKEMIN=>TURB_MODEL(KTO)%XKEMIN
+XCEDIS=>TURB_MODEL(KTO)%XCEDIS
+XCADAP=>TURB_MODEL(KTO)%XCADAP
+CTURBLEN=>TURB_MODEL(KTO)%CTURBLEN
+CTURBDIM=>TURB_MODEL(KTO)%CTURBDIM
+LTURB_FLX=>TURB_MODEL(KTO)%LTURB_FLX
+LTURB_DIAG=>TURB_MODEL(KTO)%LTURB_DIAG
+LSUBG_COND=>TURB_MODEL(KTO)%LSUBG_COND
+LSIGMAS=>TURB_MODEL(KTO)%LSIGMAS
+LSIG_CONV=>TURB_MODEL(KTO)%LSIG_CONV
+LRMC01=>TURB_MODEL(KTO)%LRMC01
+CTOM=>TURB_MODEL(KTO)%CTOM
+CSUBG_AUCV=>TURB_MODEL(KTO)%CSUBG_AUCV
+CSUBG_AUCV_RI=>TURB_MODEL(KTO)%CSUBG_AUCV_RI
+CCONDENS=>TURB_MODEL(KTO)%CCONDENS
+CLAMBDA3=>TURB_MODEL(KTO)%CLAMBDA3
+CSUBG_MF_PDF=>TURB_MODEL(KTO)%CSUBG_MF_PDF
+!XBL_DEPTH=>TURB_MODEL(KTO)%XBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
+!XSBL_DEPTH=>TURB_MODEL(KTO)%XSBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
+!XWTHVMF=>TURB_MODEL(KTO)%XWTHVMF !Done in FIELDLIST_GOTO_MODEL
+VSIGQSAT=>TURB_MODEL(KTO)%VSIGQSAT
+XDYP=>TURB_MODEL(KTO)%XDYP 
+XTHP=>TURB_MODEL(KTO)%XTHP 
+XTR=>TURB_MODEL(KTO)%XTR  
+XDISS=>TURB_MODEL(KTO)%XDISS
+XLEM=>TURB_MODEL(KTO)%XLEM
+XSSUFL_C=>TURB_MODEL(KTO)%XSSUFL_C
+XSSVFL_C=>TURB_MODEL(KTO)%XSSVFL_C
+XSSTFL_C=>TURB_MODEL(KTO)%XSSTFL_C
+XSSRFL_C=>TURB_MODEL(KTO)%XSSRFL_C
+LHGRAD=>TURB_MODEL(KTO)%LHGRAD
+XCOEFHGRADTHL=>TURB_MODEL(KTO)%XCOEFHGRADTHL
+XCOEFHGRADRM=>TURB_MODEL(KTO)%XCOEFHGRADRM
+XALTHGRAD=>TURB_MODEL(KTO)%XALTHGRAD
+XCLDTHOLD=>TURB_MODEL(KTO)%XCLDTHOLD
+
+END SUBROUTINE TURB_GOTO_MODEL
+
+END MODULE MODD_TURB_n
diff --git a/src/arome/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
similarity index 78%
rename from src/arome/turb/mode_prandtl.F90
rename to src/common/turb/mode_prandtl.F90
index d460f8e0bc4b5af58dc41c7adef903aa53ae1f67..2004edfcdf2108f3d7bc02c9ca25643aea7d1d5e 100644
--- a/src/arome/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -1,4 +1,9 @@
-!     ######spl
+!MNH_LIC Copyright 1994-2020 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
      USE PARKIND1, ONLY : JPRB
      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -6,6 +11,7 @@
 !
 !* modification 08/2010  V. Masson  smoothing of the discontinuity in functions 
 !                                   used for implicitation of exchange coefficients
+!               05/2020   V. Masson and C. Lac : bug in D_PHI3DTDZ2_O_DDTDZ
 !
 USE MODD_CTURB,      ONLY : XCTV, XCSHF, XCTD, XPHI_LIM, XCPR3, XCPR4, XCPR5
 USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
@@ -15,6 +21,539 @@ 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
+!
 SUBROUTINE SMOOTH_TURB_FUNCT(PPHI3,PF_LIM,PF)
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PPHI3   ! Phi3
diff --git a/src/common/turb/tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
similarity index 99%
rename from src/common/turb/tke_eps_sources.F90
rename to src/common/turb/mode_tke_eps_sources.F90
index f9440cfa0d7af13cd6bd8ae4feefbdcdd176c6b1..d266707a42c52b45afcf2f04ddfca6fdf3c6e506 100644
--- a/src/common/turb/tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -1,9 +1,10 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
-!-----------------------------------------------------------------
-!     ######spl
+MODULE MODE_TKE_EPS_SOURCES
+IMPLICIT NONE
+CONTAINS
       SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,  &
                     & PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,        &
                     & PTSTEP,PIMPL,PEXPL,                              &
@@ -469,3 +470,4 @@ END IF
 !
 IF (LHOOK) CALL DR_HOOK('TKE_EPS_SOURCES',1,ZHOOK_HANDLE)
 END SUBROUTINE TKE_EPS_SOURCES
+END MODULE MODE_TKE_EPS_SOURCES
diff --git a/src/arome/turb/turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
similarity index 79%
rename from src/arome/turb/turb_ver.F90
rename to src/common/turb/mode_turb_ver.F90
index 3d8b6dd0defd8f3a6c120f463490b327dfc9dcc3..32e00a0afc6051fc5a73e7608a4773a969e20e4a 100644
--- a/src/arome/turb/turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -1,9 +1,14 @@
-!     ######spl
-      SUBROUTINE TURB_VER(KKA,KKU,KKL,KRR,KRRL,KRRI,                &
-                      OCLOSE_OUT,OTURB_FLX,                         &
-                      HTURBDIM,HTOM,PIMPL,PEXPL,                    & 
-                      PTSTEP_UVW,PTSTEP_MET, PTSTEP_SV,             &
-                      HFMFILE,HLUOUT,                               &
+!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_TURB_VER
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER(KKA,KKU,KKL,KRR,KRRL,KRRI,                &
+                      OTURB_FLX,                                    &
+                      HTURBDIM,HTOM,PIMPL,PEXPL,                    &
+                      PTSTEP, TPFILE,                               &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
                       PCOSSLOPE,PSINSLOPE,                          &
                       PRHODJ,PTHVREF,                               &
@@ -15,10 +20,7 @@
                       PFWTH,PFWR,PFTH2,PFR2,PFTHR,PBL_DEPTH,        &
                       PSBL_DEPTH,PLMO,                              &
                       PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,             &
-                      PDP,PTP,PSIGS,PWTH,PWRC,PWSV)
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
+                      PDP,PTP,PSIGS,PWTH,PWRC,PWSV                  )
 !     ###############################################################
 !
 !
@@ -109,17 +111,11 @@
 !!                               field to be derivated
 !!                               _(M,UW,...) represent the localization of the 
 !!                               field	derivated
-!!                               
 !!
-!!      MXM,MXF,MYM,MYF,MZM,MZF
-!!                             :  Shuman functions (mean operators)     
-!!      DXF,DYF,DZF,DZM
-!!                             :  Shuman functions (difference operators)     
-!!                               
-!!      SUBROUTINE TRIDIAG     : to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
 !!                               of a variable located at a mass point
 !!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
 !!      FUNCTIONs ETHETA and EMOIST  :  
@@ -200,35 +196,47 @@
 !!                              advection schemes
 !!                     Feb. 2012  (Y. Seity) add possibility to run with
 !!                                 reversed vertical levels
+!!                     10/2012 (J.Escobar) Bypass PGI bug , redefine some allocatable array inplace of automatic
+!!                     08/2014 (J.Escobar) Bypass PGI memory leak bug , replace IF statement with IF THEN ENDIF
 !!      Modifications: July,    2015  (Wim de Rooy) switch for HARATU (Racmo turbulence scheme)
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! JL Redelsperger 03/2021 : add Ocean LES case
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
+USE MODD_DYN_n,          ONLY: LOCEAN
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_LES
-USE MODD_NSV, ONLY : NSV
+USE MODD_NSV,            ONLY: NSV
 !
-USE MODI_PRANDTL
+!USE MODE_PRANDTL, ONLY: PRANDTL
 USE MODI_EMOIST
 USE MODI_ETHETA
 USE MODI_GRADIENT_M
 USE MODI_GRADIENT_W
 USE MODI_TURB
-USE MODI_TURB_VER_THERMO_FLUX
-USE MODI_TURB_VER_THERMO_CORR
-USE MODI_TURB_VER_DYN_FLUX
-USE MODI_TURB_VER_SV_FLUX
-USE MODI_TURB_VER_SV_CORR
+USE MODE_TURB_VER_THERMO_FLUX, ONLY: TURB_VER_THERMO_FLUX
+USE MODE_TURB_VER_THERMO_CORR, ONLY: TURB_VER_THERMO_CORR
+USE MODE_TURB_VER_DYN_FLUX, ONLY: TURB_VER_DYN_FLUX
+USE MODE_TURB_VER_SV_FLUX, ONLY: TURB_VER_SV_FLUX
+USE MODE_TURB_VER_SV_CORR, ONLY: TURB_VER_SV_CORR
 USE MODI_LES_MEAN_SUBGRID
 USE MODI_SBL_DEPTH
+USE MODI_SECOND_MNH
 !
-USE MODE_FMWRIT
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
 USE MODE_PRANDTL
 !
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -241,21 +249,14 @@ INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -
 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.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the 
+CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
                                                       ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
+CHARACTER(len=4),       INTENT(IN)   ::  HTOM         ! type of Third Order Moment
 REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-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
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX, PDYY, PDZZ, PDZX, PDZY 
                                                       ! Metric coefficients
@@ -298,7 +299,6 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PVSLOPEM     ! wind component along the
                                      ! direction normal to the maximum slope one
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-!
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
 ! PLENGTHM PLENGTHH used in case of LHARATU
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLENGTHM     ! Turb. mixing length momentum
@@ -336,7 +336,9 @@ REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
 !
 !*       0.2  declaration of local variables
 !
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ::  &
+!JUAN BUG PGI
+!!$REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ::  &
+REAL, ALLOCATABLE, DIMENSION(:,:,:)  ::  &
        ZBETA,    & ! buoyancy coefficient
        ZSQRT_TKE,& ! sqrt(e)
        ZDTH_DZ,  & ! d(th)/dz
@@ -357,17 +359,14 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ::  &
        ZWV,      & ! (v'w')
        ZTHLP,    & ! guess of potential temperature due to vert. turbulent flux
        ZRP         ! guess of total water due to vert. turbulent flux
-REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV)  ::  &
+
+!!$REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV)  ::  &
+REAL, ALLOCATABLE, DIMENSION(:,:,:,:)  ::  &
        ZPSI_SV,  & ! Prandtl number for scalars
-       ZREDS1,   & ! 1D Redeslperger number R_sv
-       ZRED2THS, & ! 3D Redeslperger number R*2_thsv
-       ZRED2RS     ! 3D Redeslperger number R*2_rsv
+       ZREDS1,   & ! 1D Redelsperger number R_sv
+       ZRED2THS, & ! 3D Redelsperger number R*2_thsv
+       ZRED2RS     ! 3D Redelsperger number R*2_rsv
 REAL, DIMENSION(SIZE(PLM,1),SIZE(PLM,2),SIZE(PLM,3))  ::  ZLM
-INTEGER             :: IRESP        ! Return code of FM routines 
-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
 !
 LOGICAL :: GUSERV    ! flag to use water vapor
 INTEGER :: IKB,IKE   ! index value for the Beginning
@@ -375,27 +374,53 @@ INTEGER :: IKB,IKE   ! index value for the Beginning
 INTEGER :: JSV       ! loop counter on scalar variables
 REAL    :: ZTIME1
 REAL    :: ZTIME2
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+TYPE(TFIELDDATA) :: TZFIELD
+!----------------------------------------------------------------------------
+ALLOCATE (      ZBETA(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
+       ZSQRT_TKE(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)),& 
+       ZDTH_DZ(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ,&
+       ZDR_DZ(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))   ,&
+       ZRED2TH3(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ,& 
+       ZRED2R3(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ,&
+       ZRED2THR3(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)),&
+       ZBLL_O_E(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ,&
+       ZETHETA(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ,&
+       ZEMOIST(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ,&
+       ZREDTH1(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))  ,&
+       ZREDR1(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))   ,&
+       ZPHI3(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
+       ZPSI3(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
+       ZD(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))       ,&
+       ZWTHV(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
+       ZWU(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))      ,&
+       ZWV(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))      ,&
+       ZTHLP(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
+       ZRP(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))     )   
+
+ALLOCATE ( &
+ ZPSI_SV(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV),  &
+ ZREDS1(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV),   &
+ ZRED2THS(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV), &
+ ZRED2RS(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV)   )
+
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
 !             -------------
 !
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('TURB_VER',0,ZHOOK_HANDLE)
-PTP (:,:,:) = 0.
-PDP (:,:,:) = 0.
 !
 IKB=KKA+JPVEXT_TURB*KKL
 IKE=KKU-JPVEXT_TURB*KKL
-
 !
 !
 ! 3D Redelsperger numbers
 !
 !
-CALL PRANDTL(KKA,KKU,KKL,KRR,KRRI,OCLOSE_OUT,OTURB_FLX,        &
+CALL PRANDTL(KKA,KKU,KKL,KRR,KRRI,OTURB_FLX,       &
              HTURBDIM,                             &
-             HFMFILE,HLUOUT,                       &
+             TPFILE,                               &
              PDXX,PDYY,PDZZ,PDZX,PDZY,             &
              PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
              PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
@@ -404,9 +429,14 @@ CALL PRANDTL(KKA,KKU,KKL,KRR,KRRI,OCLOSE_OUT,OTURB_FLX,        &
              ZREDS1,ZRED2THS, ZRED2RS,             &
              ZBLL_O_E,                             &
              ZETHETA, ZEMOIST                      )
+!
 ! Buoyancy coefficient
 !
-ZBETA = XG/PTHVREF
+IF (LOCEAN) THEN
+  ZBETA = XG*XALPHAOC
+ELSE
+  ZBETA = XG/PTHVREF
+END IF
 !
 ! Square root of Tke
 !
@@ -422,9 +452,9 @@ IF (KRR>0) ZDR_DZ  = GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL)
 ! Denominator factor in 3rd order terms
 !
 IF (.NOT. LHARAT) THEN
-ZD(:,:,:) = (1.+ZREDTH1+ZREDR1) * (1.+0.5*(ZREDTH1+ZREDR1))
+  ZD(:,:,:) = (1.+ZREDTH1+ZREDR1) * (1.+0.5*(ZREDTH1+ZREDR1))
 ELSE
-ZD(:,:,:) = 1.
+  ZD(:,:,:) = 1.
 ENDIF
 !
 ! Phi3 and Psi3 Prandtl numbers
@@ -450,8 +480,6 @@ IF (LLES_CALL) THEN
   CALL SECOND_MNH(ZTIME2)
   XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 END IF
-
-!ENDIF
 !----------------------------------------------------------------------------
 !
 !
@@ -468,16 +496,15 @@ END IF
 !
 
 IF (LHARAT) THEN
-ZLM=PLENGTHH
+  ZLM=PLENGTHH
 ELSE
-ZLM=PLM
+  ZLM=PLM
 ENDIF
 !
   CALL  TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI,               &
-                        OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
-                        PIMPL,PEXPL,                                  &
-                        PTSTEP_MET,                                   &
-                        HFMFILE,HLUOUT,                               &
+                        OTURB_FLX,HTURBDIM,HTOM,                      &
+                        PIMPL,PEXPL,PTSTEP,                           &
+                        TPFILE,                                       &
                         PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
                         PRHODJ,PTHVREF,                               &
                         PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
@@ -491,14 +518,10 @@ ENDIF
                         MFMOIST,PBL_DEPTH,ZWTHV,                      &
                         PRTHLS,PRRS,ZTHLP,ZRP,PTP,PWTH,PWRC )
 !
-!
-!
-!  Use Lh (=Lq) from vdfexcuhl as input for turb_ver_thermo_corr
-
   CALL  TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,               &
-                        OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
+                        OTURB_FLX,HTURBDIM,HTOM,                      &
                         PIMPL,PEXPL,                                  &
-                        HFMFILE,HLUOUT,                               &
+                        TPFILE,                                       &
                         PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,           &
                         PRHODJ,PTHVREF,                               &
                         PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
@@ -525,15 +548,12 @@ ENDIF
 !             -----------------------------------------------
 !
 !
-IF (LHARAT) THEN
-ZLM=PLENGTHM
-ENDIF
-
+IF (LHARAT) ZLM=PLENGTHM
+!
 CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
-                      OCLOSE_OUT,OTURB_FLX,KRR,                     &
-                      HTURBDIM,PIMPL,PEXPL,                         &
-                      PTSTEP_UVW,                                   &
-                      HFMFILE,HLUOUT,                               &
+                      OTURB_FLX,KRR,                                &
+                      HTURBDIM,PIMPL,PEXPL,PTSTEP,                  &
+                      TPFILE,                                       &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
                       PCOSSLOPE,PSINSLOPE,                          &
                       PRHODJ,                                       &
@@ -541,7 +561,7 @@ CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
                       PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
                       PTKEM,ZLM,MFMOIST,ZWU,ZWV,                    &
                       PRUS,PRVS,PRWS,                               &
-                      PDP,PTP                               )
+                      PDP,PTP                                       )
 !
 !----------------------------------------------------------------------------
 !
@@ -549,16 +569,13 @@ CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
 !*       8.   SOURCES OF PASSIVE SCALAR VARIABLES
 !             -----------------------------------
 !
-IF (LHARAT) THEN
-ZLM=PLENGTHH
-ENDIF
-
+IF (LHARAT) ZLM=PLENGTHH
+!
 IF (SIZE(PSVM,4)>0)                                                 &
 CALL  TURB_VER_SV_FLUX(KKA,KKU,KKL,                                 &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,                &
-                      PIMPL,PEXPL,                                  &
-                      PTSTEP_SV,                                    &
-                      HFMFILE,HLUOUT,                               &
+                      OTURB_FLX,HTURBDIM,                           &
+                      PIMPL,PEXPL,PTSTEP,                           &
+                      TPFILE,                                       &
                       PDZZ,PDIRCOSZW,                               &
                       PRHODJ,PWM,                                   &
                       PSFSVM,PSFSVP,                                &
@@ -568,7 +585,7 @@ CALL  TURB_VER_SV_FLUX(KKA,KKU,KKL,                                 &
 !
 !
 IF (SIZE(PSVM,4)>0 .AND. LLES_CALL)                                 &
-CALL  TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,                               &
+CALL  TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,                   &
                       PDZZ,                                         &
                       PTHLM,PRM,PTHVREF,                            &
                       PLOCPEXNM,PATHETA,PAMOIST,PSRCM,ZPHI3,ZPSI3,  &
@@ -590,40 +607,57 @@ IF (SIZE(PSBL_DEPTH)>0) CALL SBL_DEPTH(IKB,IKE,PZZ,ZWU,ZWV,ZWTHV,PLMO,PSBL_DEPTH
 !             ------
 !
 !
-IF ( OTURB_FLX .AND. OCLOSE_OUT .AND. .NOT. LHARAT) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. .NOT. LHARAT) THEN
 !
 ! stores the Turbulent Prandtl number
 ! 
-  YRECFM  ='PHI3'
-  YCOMMENT='X_Y_Z_PHI3 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-  CALL  FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZPHI3,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'PHI3'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'PHI3'
+  TZFIELD%CUNITS     = '1'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'Turbulent Prandtl number'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZPHI3)
 !
 ! stores the Turbulent Schmidt number
 ! 
-  YRECFM  ='PSI3'
-  YCOMMENT='X_Y_Z_PSI3 (0)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)
-!
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZPSI3,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'PSI3'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'PSI3'
+  TZFIELD%CUNITS     = '1'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'Turbulent Schmidt number'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZPSI3)
 !
 !
 ! stores the Turbulent Schmidt number for the scalar variables
 ! 
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CUNITS     = '1'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
   DO JSV=1,NSV
-    WRITE(YRECFM, '("PSI_SV_",I3.3)') JSV
-    YCOMMENT='X_Y_Z_'//YRECFM//' (0)'
-    IGRID   = 4
-    ILENCH=LEN(YCOMMENT)
-    CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZPSI_SV(:,:,:,JSV),   &
-                IGRID,ILENCH,YCOMMENT,IRESP)
+    WRITE(TZFIELD%CMNHNAME, '("PSI_SV_",I3.3)') JSV
+    TZFIELD%CLONGNAME  = TRIM(TZFIELD%CMNHNAME)
+    TZFIELD%CCOMMENT   = 'X_Y_Z_'//TRIM(TZFIELD%CMNHNAME)
+    CALL IO_Field_write(TPFILE,TZFIELD,ZPSI_SV(:,:,:,JSV))
   END DO
-
+!
 END IF
 !
 !
 !----------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('TURB_VER',1,ZHOOK_HANDLE)
-END SUBROUTINE TURB_VER                                                                                               
+END SUBROUTINE TURB_VER
+END MODULE MODE_TURB_VER 
diff --git a/src/arome/turb/turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
similarity index 91%
rename from src/arome/turb/turb_ver_dyn_flux.F90
rename to src/common/turb/mode_turb_ver_dyn_flux.F90
index 6c3c5f5c7504fcf51cf7bb70a23eddc40539b42b..4473628de977d93638286ef291b179b7dc68e47a 100644
--- a/src/arome/turb/turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -1,9 +1,15 @@
-!     ######spl
-      SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL,                     &
-                      OCLOSE_OUT,OTURB_FLX,KRR,                     &
+!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_TURB_VER_DYN_FLUX
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL,                     &
+                      OTURB_FLX,KRR,                                &
                       HTURBDIM,PIMPL,PEXPL,                         &
                       PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
+                      TPFILE,                                       &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
                       PCOSSLOPE,PSINSLOPE,                          &
                       PRHODJ,                                       &
@@ -12,10 +18,6 @@
                       PTKEM,PLM,MFMOIST,PWU,PWV,                    &
                       PRUS,PRVS,PRWS,                               &
                       PDP,PTP                                       )
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
-
 !     ###############################################################
 !
 !
@@ -113,10 +115,10 @@
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
 !!                               of a variable located at a mass point
 !!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
 !!      FUNCTIONs ETHETA and EMOIST  :  
@@ -195,31 +197,46 @@
 !!                                              change of YCOMMENT
 !!      2012-02 Y. Seity,  add possibility to run with reversed vertical levels
 !!      Modifications  July 2015 (Wim de Rooy) LHARATU switch
-
+!!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!      Q. Rodier      17/01/2019 : cleaning : remove cyclic conditions on DP and ZA
+!! JL Redelsperger 03/2021 : Add Ocean  & O-A Autocoupling LES Cases
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CONF
 USE MODD_CST
 USE MODD_CTURB
-USE MODD_PARAMETERS
+USE MODD_DYN_n,          ONLY: LOCEAN
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_LES
 USE MODD_NSV
+USE MODD_OCEANH
+USE MODD_PARAMETERS
+USE MODD_REF, ONLY : LCOUPLES
+USE MODD_TURB_n
 !
 !
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
+USE MODI_SECOND_MNH
 USE MODI_SHUMAN , ONLY: MZM, MZF, MXM, MXF, MYM, MYF,&
                       & DZM, DXF, DXM, DYF, DYM
 USE MODI_TRIDIAG 
 USE MODI_TRIDIAG_WIND 
-USE MODE_FMWRIT
 USE MODI_LES_MEAN_SUBGRID
 !
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+USE MODE_ll
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -229,19 +246,14 @@ IMPLICIT NONE
 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
-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
 INTEGER,                INTENT(IN)   ::  KRR          ! number of moist var.
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the 
+CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
                                                       ! turbulence scheme
 REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
 REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file 
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDXX, PDYY, PDZZ, PDZX, PDZY 
                                                       ! Metric coefficients
@@ -281,8 +293,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PWV          ! momentum flux v'w'
 REAL, DIMENSION(:,:,:), INTENT(INOUT)   ::  PRUS, PRVS, PRWS
                             ! cumulated sources for the prognostic variables
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT)::  PDP,PTP   ! Dynamic and thermal
-                                                   ! TKE production terms
+REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PDP          ! Dynamic TKE production term
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTP          ! Thermal TKE production term
 !
 !
 !
@@ -323,6 +335,7 @@ REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
 INTEGER             :: IIU,IJU      ! size of array in x,y,z directions
 !
 REAL :: ZTIME1, ZTIME2, ZCMFS
+TYPE(TFIELDDATA) :: TZFIELD
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
@@ -444,13 +457,19 @@ ZFLXZ(:,:,IKB:IKB)   =   MXM(PDZZ(:,:,IKB:IKB))  *                &
 ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB) 
 
 !
-IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
   ! stores the U wind component vertical flux
-  YRECFM  ='UW_VFLX'
-  YCOMMENT='X_Y_Z_UW_VFLX (M**2/S**2)'
-  IGRID   = 4  
-  ILENCH=LEN(YCOMMENT) 
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'UW_VFLX'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'UW_VFLX'
+  TZFIELD%CUNITS     = 'm2 s-2'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'U wind component vertical flux'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
 END IF
 !
 ! first part of total momentum flux
@@ -613,13 +632,19 @@ ZFLXZ(:,:,IKB:IKB)   =   MYM(PDZZ(:,:,IKB:IKB))  *                       &
 !
 ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
 !
-IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
   ! stores the V wind component vertical flux
-  YRECFM  ='VW_VFLX'
-  YCOMMENT='X_Y_Z_VW_VFLX (M**2/S**2)'
-  IGRID   = 4
-  ILENCH=LEN(YCOMMENT)  
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'VW_VFLX'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'VW_VFLX'
+  TZFIELD%CUNITS     = 'm2 s-2'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'V wind component vertical flux'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
 END IF
 !
 ! second part of total momentum flux
@@ -727,20 +752,27 @@ END IF
 !*       7.   DIAGNOSTIC COMPUTATION OF THE 1D <W W> VARIANCE
 !             -----------------------------------------------
 !
-IF ( OTURB_FLX .AND. OCLOSE_OUT .AND. HTURBDIM == '1DIM') THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. HTURBDIM == '1DIM') THEN
   ZFLXZ(:,:,:)= (2./3.) * PTKEM(:,:,:)                     &
      -ZCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)
   ! to be tested &
   !   +XCMFB*(4./3.)*PLM(:,:,:)/SQRT(PTKEM(:,:,:))*PTP(:,:,:) 
   ! stores the W variance
-  YRECFM  ='W_VVAR'
-  YCOMMENT='X_Y_Z_W_VVAR (M**2/S**2)'
-  IGRID   = 1  
-  ILENCH=LEN(YCOMMENT) 
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'W_VVAR'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'W_VVAR'
+  TZFIELD%CUNITS     = 'm2 s-2'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'X_Y_Z_W_VVAR'
+  TZFIELD%NGRID      = 1
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
 END IF
 !
 !----------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_DYN_FLUX
+END MODULE MODE_TURB_VER_DYN_FLUX
diff --git a/src/arome/turb/turb_ver_sv_corr.F90 b/src/common/turb/mode_turb_ver_sv_corr.F90
similarity index 91%
rename from src/arome/turb/turb_ver_sv_corr.F90
rename to src/common/turb/mode_turb_ver_sv_corr.F90
index 25744253b913e3b7ea2203fb5e77e26b6de466ce..f140c0792e47cc3315f874e56bd9f2241598d69e 100644
--- a/src/arome/turb/turb_ver_sv_corr.F90
+++ b/src/common/turb/mode_turb_ver_sv_corr.F90
@@ -1,12 +1,16 @@
-!     ######spl
-      SUBROUTINE TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,        &
+!MNH_LIC Copyright 2002-2020 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_TURB_VER_SV_CORR
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,        &
                       PDZZ,                                         &
                       PTHLM,PRM,PTHVREF,                            &
                       PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PPHI3,PPSI3,  &
                       PWM,PSVM,                                     &
                       PTKEM,PLM,PLEPS,PPSI_SV                       )
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !     ###############################################################
 !
 !
@@ -48,12 +52,16 @@
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_PARAMETERS
 USE MODD_LES
 USE MODD_CONF
 USE MODD_NSV, ONLY : NSV,NSV_LGBEG,NSV_LGEND
+USE MODD_BLOWSNOW
 !
 !
 USE MODI_GRADIENT_U
@@ -65,6 +73,8 @@ USE MODI_EMOIST
 USE MODI_ETHETA
 USE MODI_LES_MEAN_SUBGRID
 !
+USE MODI_SECOND_MNH
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -104,6 +114,9 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PPSI_SV      ! Inv.Turb.Sch.for scalars
 !
 REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3))  ::  &
        ZA, ZFLXZ
+!
+REAL :: ZCSV          !constant for the scalar flux
+!
 INTEGER             :: JSV          ! loop counters
 !
 REAL :: ZTIME1, ZTIME2
@@ -162,3 +175,4 @@ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 !
 IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_SV_CORR
+END MODULE MODE_TURB_VER_SV_CORR
diff --git a/src/arome/turb/turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
similarity index 88%
rename from src/arome/turb/turb_ver_sv_flux.F90
rename to src/common/turb/mode_turb_ver_sv_flux.F90
index d19712dc51381e9a0f56ace65d0f65e95bdb98dc..9ef220cf8e7ca6883f9e689bcfeaa8cdfd3041bf 100644
--- a/src/arome/turb/turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -1,9 +1,15 @@
-!     ######spl
-      SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL,                      &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,                &
+!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_TURB_VER_SV_FLUX
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL,                      &
+                      OTURB_FLX,HTURBDIM,                           &
                       PIMPL,PEXPL,                                  &
                       PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
+                      TPFILE,                                       &
                       PDZZ,PDIRCOSZW,                               &
                       PRHODJ,PWM,                                   &
                       PSFSVM,PSFSVP,                                &
@@ -11,10 +17,6 @@
                       PTKEM,PLM,MFMOIST,PPSI_SV,                    &
                       PRSVS,PWSV                                    )
 !
-
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
 !
 !
 !!****  *TURB_VER_SV_FLUX* -compute the source terms due to the vertical turbulent
@@ -111,10 +113,10 @@
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
 !!                               of a variable located at a mass point
 !!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
 !!      FUNCTIONs ETHETA and EMOIST  :  
@@ -195,18 +197,29 @@
 !!                                              change of YCOMMENT
 !!                     Feb 2012(Y. Seity) add possibility to run with reversed 
 !!                                              vertical levels
-!!      Modifications: July 2015 (Wim de Rooy) LHARATU switch
+!!      Modifications: July 2015 (Wim de Rooy) LHARAT switch
+!!                     Feb 2017(M. Leriche) add initialisation of ZSOURCE
+!!                                   to avoid unknwon values outside physical domain
+!!                                   and avoid negative values in sv tendencies
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_LES
 USE MODD_CONF
-USE MODD_NSV, ONLY : NSV_LGBEG,NSV_LGEND
+USE MODD_NSV,            ONLY: XSVMIN, NSV_LGBEG, NSV_LGEND
+USE MODD_BLOWSNOW
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 !
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
@@ -217,9 +230,10 @@ USE MODI_TRIDIAG
 USE MODI_TRIDIAG_WIND 
 USE MODI_EMOIST
 USE MODI_ETHETA
-USE MODE_FMWRIT
 USE MODI_LES_MEAN_SUBGRID
 !
+USE MODI_SECOND_MNH
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -228,18 +242,13 @@ IMPLICIT NONE
 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
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
+CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
                                                       ! turbulence scheme
 REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
 REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file 
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDZZ
                                                       ! Metric coefficients
@@ -297,6 +306,9 @@ CHARACTER (LEN=16)  :: YRECFM       ! Name of the desired field in LFIFM file
 REAL :: ZTIME1, ZTIME2
 
 REAL :: ZCSVP = 4.0  ! constant for scalar flux presso-correlation (RS81)
+REAL :: ZCSV          !constant for the scalar flux
+!
+TYPE(TFIELDDATA)  :: TZFIELD
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
@@ -359,14 +371,14 @@ ENDIF
   ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
   ZSOURCE(:,:,IKE) = 0.
 !
-! Obtention of the splitted JSV scalar variable at t+ deltat  
+! Obtention of the split JSV scalar variable at t+ deltat  
   CALL TRIDIAG(KKA,KKU,KKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
 !
 !  Compute the equivalent tendency for the JSV scalar variable
   PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV)+    &
                     PRHODJ(:,:,:)*(ZRES(:,:,:)-PSVM(:,:,:,JSV))/PTSTEP
 !
-  IF ( (OTURB_FLX .AND. OCLOSE_OUT) .OR. LLES_CALL ) THEN
+  IF ( (OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL ) THEN
     ! Diagnostic of the cartesian vertical flux
     !
     ZFLXZ(:,:,:) = -XCHF * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), KKA, KKU, KKL) / PDZZ * &
@@ -393,17 +405,25 @@ ENDIF
     PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-KKL,JSV)
  END IF
   !
-  IF (OTURB_FLX .AND. OCLOSE_OUT) THEN
+  IF (OTURB_FLX .AND. TPFILE%LOPENED) THEN
     ! stores the JSVth vertical flux
-    WRITE(YRECFM,'("WSV_FLX_",I3.3)') JSV 
-    YCOMMENT='X_Y_Z_'//YRECFM//' (SVUNIT*M/S)'
-    IGRID   = 4  
-    ILENCH=LEN(YCOMMENT)
-    CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+    WRITE(TZFIELD%CMNHNAME,'("WSV_FLX_",I3.3)') JSV
+    TZFIELD%CSTDNAME   = ''
+    TZFIELD%CLONGNAME  = TRIM(TZFIELD%CMNHNAME)
+    !PW: TODO: use the correct units of the JSV variable (and multiply it by m s-1)
+    TZFIELD%CUNITS     = 'SVUNIT m s-1'
+    TZFIELD%CDIR       = 'XY'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_'//TRIM(TZFIELD%CMNHNAME)
+    TZFIELD%NGRID      = 4
+    TZFIELD%NTYPE      = TYPEREAL
+    TZFIELD%NDIMS      = 3
+    TZFIELD%LTIMEDEP   = .TRUE.
+    !
+    CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
   END IF
   !
   ! Storage in the LES configuration
-  
+  !
   IF (LLES_CALL) THEN
     CALL SECOND_MNH(ZTIME1)
     CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WSv(:,:,:,JSV) )
@@ -423,3 +443,4 @@ END DO   ! end of scalar loop
 !
 IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_SV_FLUX
+END MODULE MODE_TURB_VER_SV_FLUX
diff --git a/src/arome/turb/turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
similarity index 93%
rename from src/arome/turb/turb_ver_thermo_corr.F90
rename to src/common/turb/mode_turb_ver_thermo_corr.F90
index 1e343bb2705068a5d5dc6dd48f38c7a405b577ec..3d420f3936d33e208d20b38c8e15974356c54d1d 100644
--- a/src/arome/turb/turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -1,8 +1,14 @@
-!     ######spl
-      SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
+!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_TURB_VER_THERMO_CORR
+IMPLICIT NONE
+CONTAINS      
+SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
+                      OTURB_FLX,HTURBDIM,HTOM,                      &
                       PIMPL,PEXPL,                                  &
-                      HFMFILE,HLUOUT,                               &
+                      TPFILE,                                       &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,           &
                       PRHODJ,PTHVREF,                               &
                       PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
@@ -13,11 +19,7 @@
                       PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
                       PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
                       PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
-                      PTHLP,PRP,MFMOIST,PSIGS                  )
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
-
+                      PTHLP,PRP,MFMOIST,PSIGS                       )
 !     ###############################################################
 !
 !
@@ -115,10 +117,10 @@
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
 !!                               of a variable located at a mass point
 !!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
 !!      FUNCTIONs ETHETA and EMOIST  :  
@@ -147,7 +149,7 @@
 !!
 !!      Module MODD_PARAMETERS
 !!
-!!           JPVEXT     : number of vertical external points
+!!           JPVEXT_TURB     : number of vertical external points
 !!           JPHEXT     : number of horizontal external points
 !!
 !!
@@ -199,14 +201,19 @@
 !!                                              change of YCOMMENT
 !!                     2012-02 (Y. Seity) add possibility to run with reversed 
 !!                                              vertical levels
-!!      Modifications  July 2015 (Wim de Rooy) LHARATU switch
+!!      Modifications  July 2015 (Wim de Rooy) LHARAT switch
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 USE MODD_CST
 USE MODD_CTURB
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_CONF
 USE MODD_LES
@@ -217,13 +224,14 @@ USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
 USE MODI_TRIDIAG 
-USE MODE_FMWRIT
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_PRANDTL
 USE MODI_TRIDIAG_THERMO
 !
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
 USE MODE_PRANDTL
 !
+USE MODI_SECOND_MNH
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -236,18 +244,13 @@ INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -
 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.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
+CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
                                                       ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
+CHARACTER(len=4),       INTENT(IN)   ::  HTOM         ! type of Third Order Moment
 REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file 
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDZZ, PDXX, PDYY, PDZX, PDZY
                                                       ! Metric coefficients
@@ -349,7 +352,7 @@ LOGICAL :: GFWTH    ! flag to use w'2th'
 LOGICAL :: GFR2     ! flag to use w'r'2
 LOGICAL :: GFWR     ! flag to use w'2r'
 LOGICAL :: GFTHR    ! flag to use w'th'r'
-
+TYPE(TFIELDDATA) :: TZFIELD
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
@@ -513,12 +516,18 @@ ENDIF
   !
   !
   ! stores <THl THl>  
-  IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
-    YRECFM  ='THL_VVAR'
-    YCOMMENT='X_Y_Z_THL_VVAR (KELVIN**2)'
-    IGRID   = 1
-    ILENCH=LEN(YCOMMENT)
-    CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+    TZFIELD%CMNHNAME   = 'THL_VVAR'
+    TZFIELD%CSTDNAME   = ''
+    TZFIELD%CLONGNAME  = 'THL_VVAR'
+    TZFIELD%CUNITS     = 'K2'
+    TZFIELD%CDIR       = 'XY'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_THL_VVAR'
+    TZFIELD%NGRID      = 1
+    TZFIELD%NTYPE      = TYPEREAL
+    TZFIELD%NDIMS      = 3
+    TZFIELD%LTIMEDEP   = .TRUE.
+    CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
   END IF
 !
 ! and we store in LES configuration
@@ -674,12 +683,18 @@ ENDIF
                      2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLXZ(:,:,:)
     END IF
     ! stores <THl Rnp>   
-    IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
-      YRECFM  ='THLRCONS_VCOR'
-      YCOMMENT='X_Y_Z_THLRCONS_VCOR (KELVIN*KG/KG)'
-      IGRID   = 1
-      ILENCH=LEN(YCOMMENT)
-      CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+    IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+      TZFIELD%CMNHNAME   = 'THLRCONS_VCOR'
+      TZFIELD%CSTDNAME   = ''
+      TZFIELD%CLONGNAME  = 'THLRCONS_VCOR'
+      TZFIELD%CUNITS     = 'K kg kg-1'
+      TZFIELD%CDIR       = 'XY'
+      TZFIELD%CCOMMENT   = 'X_Y_Z_THLRCONS_VCOR'
+      TZFIELD%NGRID      = 1
+      TZFIELD%NTYPE      = TYPEREAL
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .TRUE.
+      CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
     END IF
 !
 ! and we store in LES configuration
@@ -801,12 +816,18 @@ ENDIF
       PSIGS(:,:,:) = PSIGS(:,:,:) + PAMOIST(:,:,:) **2 * ZFLXZ(:,:,:)
     END IF
     ! stores <Rnp Rnp>    
-    IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
-      YRECFM  ='RTOT_VVAR'
-      YCOMMENT='X_Y_Z_RTOT_VVAR (KG/KG **2)'
-      IGRID   = 1
-      ILENCH=LEN(YCOMMENT)
-      CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+    IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+      TZFIELD%CMNHNAME   = 'RTOT_VVAR'
+      TZFIELD%CSTDNAME   = ''
+      TZFIELD%CLONGNAME  = 'RTOT_VVAR'
+      TZFIELD%CUNITS     = 'kg2 kg-2'
+      TZFIELD%CDIR       = 'XY'
+      TZFIELD%CCOMMENT   = 'X_Y_Z_RTOT_VVAR'
+      TZFIELD%NGRID      = 1
+      TZFIELD%NTYPE      = TYPEREAL
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .TRUE.
+      CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
     END IF
     !
     ! and we store in LES configuration
@@ -842,3 +863,4 @@ ENDIF
 !----------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_CORR',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_THERMO_CORR
+END MODULE MODE_TURB_VER_THERMO_CORR
diff --git a/src/arome/turb/turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
similarity index 89%
rename from src/arome/turb/turb_ver_thermo_flux.F90
rename to src/common/turb/mode_turb_ver_thermo_flux.F90
index d377e2c9d3de1ec857dfe3167d1e89c218d3c6c6..7a79162a45544c8ea8ebb7eab1e52629cad69006 100644
--- a/src/arome/turb/turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -1,9 +1,16 @@
-!     ######spl
-      SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
-                      OCLOSE_OUT,OTURB_FLX,HTURBDIM,HTOM,           &
+!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_TURB_VER_THERMO_FLUX
+IMPLICIT NONE
+CONTAINS
+      
+SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
+                      OTURB_FLX,HTURBDIM,HTOM,                      &
                       PIMPL,PEXPL,                                  &
                       PTSTEP,                                       &
-                      HFMFILE,HLUOUT,                               &
+                      TPFILE,                                       &
                       PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,       &
                       PRHODJ,PTHVREF,                               &
                       PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
@@ -14,10 +21,7 @@
                       PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
                       PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
                       PFWTH,PFWR,PFTH2,PFR2,PFTHR,MFMOIST,PBL_DEPTH,&
-                      PWTHV,PRTHLS,PRRS,PTHLP,PRP,PTP,PWTH,PWRC)
-      USE PARKIND1, ONLY : JPRB
-      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-      USE MODD_CTURB, ONLY : LHARAT
+                      PWTHV,PRTHLS,PRRS,PTHLP,PRP,PTP,PWTH,PWRC     )
 !     ###############################################################
 !
 !
@@ -115,10 +119,10 @@
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
 !!                               of a variable located at a mass point
 !!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
+!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
 !!      FUNCTIONs ETHETA and EMOIST  :  
@@ -206,16 +210,39 @@
 !!                     2012-02 (Y. Seity) add possibility to run with reversed
 !!                                             vertical levels
 !!      Modifications  July 2015 (Wim de Rooy) LHARAT switch
+!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!                     2021 (D. Ricard) last version of HGRAD turbulence scheme
+!!                                 Leronard terms instead of Reynolds terms
+!!                                 applied to vertical fluxes of r_np and Thl
+!!                                 for implicit version of turbulence scheme
+!!                                 corrections and cleaning
+!!                     June 2020 (B. Vie) Patch preventing negative rc and ri in 2.3 and 3.3
+!! JL Redelsperger  : 03/2021: Ocean and Autocoupling O-A LES Cases
+!!                             Sfc flux shape for LDEEPOC Case
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
+USE MODD_GRID_n,         ONLY: XZS, XXHAT, XYHAT
+USE MODD_IO,             ONLY: TFILEDATA
+USE MODD_METRICS_n,      ONLY: XDXX, XDYY, XDZX, XDZY, XDZZ
 USE MODD_PARAMETERS
+USE MODD_TURB_n,         ONLY: LHGRAD, XCOEFHGRADTHL, XCOEFHGRADRM, XALTHGRAD, XCLDTHOLD
 USE MODD_CONF
 USE MODD_LES
+USE MODD_DIM_n
+USE MODD_DYN_n,          ONLY: LOCEAN
+USE MODD_OCEANH
+USE MODD_REF,            ONLY: LCOUPLES
+USE MODD_TURB_n
+USE MODD_FRC
 !
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
@@ -223,14 +250,17 @@ USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM
 USE MODI_TRIDIAG 
-USE MODE_FMWRIT
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_PRANDTL
 USE MODI_TRIDIAG_THERMO
 USE MODI_TM06_H
 !
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_PRANDTL
 !
+USE MODI_SECOND_MNH
+USE MODE_ll
+USE MODE_GATHER_ll
+!
 IMPLICIT NONE
 !
 !*      0.1  declarations of arguments
@@ -243,19 +273,14 @@ INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -
 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.
-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
-CHARACTER*4,            INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
+CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
                                                       ! turbulence scheme
-CHARACTER*4,            INTENT(IN)   ::  HTOM         ! type of Third Order Moment
+CHARACTER(len=4),       INTENT(IN)   ::  HTOM         ! type of Third Order Moment
 REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
 REAL,                   INTENT(IN)   ::  PTSTEP       ! Double Time Step
-CHARACTER(LEN=*),       INTENT(IN)   ::  HFMFILE      ! Name of the output
-                                                      ! FM-file 
-CHARACTER(LEN=*),       INTENT(IN)   ::  HLUOUT       ! Output-listing name for
-                                                      ! model n
+TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDZZ, PDXX, PDYY, PDZX, PDZY
                                                       ! Metric coefficients
@@ -321,7 +346,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS       ! cumulated source for rt
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTHLP      ! guess of thl at t+ deltat
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRP        ! guess of r at t+ deltat
 !
-REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PTP       ! Dynamic and thermal
+REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTP       ! Dynamic and thermal
                                                      ! TKE production terms
 !
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PWTH       ! heat flux
@@ -359,6 +384,7 @@ LOGICAL :: GFWTH    ! flag to use w'2th'
 LOGICAL :: GFR2     ! flag to use w'r'2
 LOGICAL :: GFWR     ! flag to use w'2r'
 LOGICAL :: GFTHR    ! flag to use w'th'r'
+TYPE(TFIELDDATA) :: TZFIELD
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
@@ -503,13 +529,19 @@ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
   PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
   PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
 
-IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
   ! stores the conservative potential temperature vertical flux
-  YRECFM  ='THW_FLX'
-  YCOMMENT='X_Y_Z_THW_FLX (K*M/S)'
-  IGRID   = 4  
-  ILENCH=LEN(YCOMMENT) 
-  CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  TZFIELD%CMNHNAME   = 'THW_FLX'
+  TZFIELD%CSTDNAME   = ''
+  TZFIELD%CLONGNAME  = 'THW_FLX'
+  TZFIELD%CUNITS     = 'K m s-1'
+  TZFIELD%CDIR       = 'XY'
+  TZFIELD%CCOMMENT   = 'Conservative potential temperature vertical flux'
+  TZFIELD%NGRID      = 4
+  TZFIELD%NTYPE      = TYPEREAL
+  TZFIELD%NDIMS      = 3
+  TZFIELD%LTIMEDEP   = .TRUE.
+  CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
 END IF
 !
 ! Contribution of the conservative temperature flux to the buoyancy flux
@@ -680,13 +712,19 @@ IF (KRR /= 0) THEN
   PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
   !
   !
-  IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
     ! stores the conservative mixing ratio vertical flux
-    YRECFM  ='RCONSW_FLX'
-    YCOMMENT='X_Y_Z_RCONSW_FLX (KG*M/S/KG)'
-    IGRID   = 4  
-    ILENCH=LEN(YCOMMENT) 
-    CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+    TZFIELD%CMNHNAME   = 'RCONSW_FLX'
+    TZFIELD%CSTDNAME   = ''
+    TZFIELD%CLONGNAME  = 'RCONSW_FLX'
+    TZFIELD%CUNITS     = 'kg m s-1 kg-1'
+    TZFIELD%CDIR       = 'XY'
+    TZFIELD%CCOMMENT   = 'Conservative mixing ratio vertical flux'
+    TZFIELD%NGRID      = 4
+    TZFIELD%NTYPE      = TYPEREAL
+    TZFIELD%NDIMS      = 3
+    TZFIELD%LTIMEDEP   = .TRUE.
+    CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
   END IF
   !
   ! Contribution of the conservative water flux to the Buoyancy flux
@@ -743,7 +781,7 @@ END IF
 !
 !*       4.1  <w Rc>    
 !
-IF ( ((OTURB_FLX .AND. OCLOSE_OUT) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
+IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
 !  
 ! recover the Conservative potential temperature flux : 
 ! With LHARAT is true tke and length scales at half levels
@@ -764,12 +802,18 @@ IF ( ((OTURB_FLX .AND. OCLOSE_OUT) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
   ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB) 
   !                 
   ! store the liquid water mixing ratio vertical flux
-  IF ( OTURB_FLX .AND. OCLOSE_OUT ) THEN
-    YRECFM  ='RCW_FLX'
-    YCOMMENT='X_Y_Z_RCW_FLX (KG*M/S/KG)'
-    IGRID   = 4  
-    ILENCH=LEN(YCOMMENT) 
-    CALL FMWRIT(HFMFILE,YRECFM,HLUOUT,'XY',ZFLXZ,IGRID,ILENCH,YCOMMENT,IRESP)
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+    TZFIELD%CMNHNAME   = 'RCW_FLX'
+    TZFIELD%CSTDNAME   = ''
+    TZFIELD%CLONGNAME  = 'RCW_FLX'
+    TZFIELD%CUNITS     = 'kg m s-1 kg-1'
+    TZFIELD%CDIR       = 'XY'
+    TZFIELD%CCOMMENT   = 'Liquid water mixing ratio vertical flux'
+    TZFIELD%NGRID      = 4
+    TZFIELD%NTYPE      = TYPEREAL
+    TZFIELD%NDIMS      = 3
+    TZFIELD%LTIMEDEP   = .TRUE.
+    CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
   END IF
   !  
 ! and we store in LES configuration this subgrid flux <w'rc'>
@@ -786,3 +830,4 @@ END IF !end of <w Rc>
 !----------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_THERMO_FLUX
+END MODULE MODE_TURB_VER_THERMO_FLUX
diff --git a/src/common/turb/prandtl.F90 b/src/common/turb/prandtl.F90
new file mode 100644
index 0000000000000000000000000000000000000000..458899f3908ef3d12ae9a65efd25f416b8245f46
--- /dev/null
+++ b/src/common/turb/prandtl.F90
@@ -0,0 +1,540 @@
+!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 04af09a38ec1e638eddd91f6179647b5cdcf3dae..66b6b86f4804a50f8e00e76ba3c93a76b4b98989 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -230,11 +230,11 @@ USE MODD_LES
 USE MODD_NSV
 !
 USE MODI_BL89
-USE MODI_TURB_VER
+USE MODE_TURB_VER, ONLY : TURB_VER
 !!MODIF AROME
 !USE MODI_ROTATE_WIND
 !USE MODI_TURB_HOR_SPLT 
-USE MODI_TKE_EPS_SOURCES
+USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
 USE MODI_SHUMAN, ONLY : MZF, MXF, MYF
 USE MODI_GRADIENT_M
 USE MODI_BUDGET_DDH
@@ -808,10 +808,9 @@ ZFTHR(:,:,:IKTB) = 0.
 !          -----------------
 !
 CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI,               &
-          OCLOSE_OUT,OTURB_FLX,                          &
+          OTURB_FLX,                                     &
           HTURBDIM,HTOM,PIMPL,ZEXPL,                     &
-          PTSTEP_UVW, PTSTEP_MET, PTSTEP_SV,             &
-          HFMFILE,HLUOUT,                                &
+          PTSTEP_MET,TPFILE,                                 &
           PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,        &
           PCOSSLOPE,PSINSLOPE,                           &
           PRHODJ,PTHVREF,                                &
diff --git a/src/mesonh/turb/prandtl.f90 b/src/mesonh/turb/prandtl.f90
index fbfe0a7621714cebb151faee288c99338a4a555b..9b8455f87638856b4ee123aa972f1223ba0a4754 100644
--- a/src/mesonh/turb/prandtl.f90
+++ b/src/mesonh/turb/prandtl.f90
@@ -185,6 +185,7 @@ END MODULE MODI_PRANDTL
 !!                                              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 
@@ -221,7 +222,7 @@ 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.
+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
@@ -535,7 +536,7 @@ END IF ! end of HTURBDIM if-block
 !*          6. SAVES THE REDELSPERGER NUMBERS
 !              ------------------------------
 !
-IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN
+IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
   !
   ! stores the RED_TH1
   TZFIELD%CMNHNAME   = 'RED_TH1'
@@ -606,4 +607,5 @@ END IF
 !
 !---------------------------------------------------------------------------
 !
+IF (LHOOK) CALL DR_HOOK('PRANDTL',1,ZHOOK_HANDLE)
 END SUBROUTINE PRANDTL
diff --git a/src/mesonh/turb/turb_ver.f90 b/src/mesonh/turb/turb_ver.f90
index 4117d8191eb9def704b654e606eba90307fe65ec..420617b508017a6f3f6d311b5d4bddcf9ae11750 100644
--- a/src/mesonh/turb/turb_ver.f90
+++ b/src/mesonh/turb/turb_ver.f90
@@ -3,12 +3,7 @@
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!    #################### 
-     MODULE MODI_TURB_VER 
-!    ####################
-!
-INTERFACE 
-!
+!     ######spl
       SUBROUTINE TURB_VER(KKA,KKU,KKL,KRR,KRRL,KRRI,                &
                       OTURB_FLX,                                    &
                       HTURBDIM,HTOM,PIMPL,PEXPL,                    &
@@ -19,7 +14,7 @@ INTERFACE
                       PSFTHM,PSFRM,PSFSVM,PSFTHP,PSFRP,PSFSVP,      &
                       PCDUEFF,PTAU11M,PTAU12M,PTAU33M,              &
                       PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
-                      PTKEM,PLM,PLEPS,                              &
+                      PTKEM,PLM,PLENGTHM,PLENGTHH,PLEPS,MFMOIST,    &
                       PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE,    &
                       PFWTH,PFWR,PFTH2,PFR2,PFTHR,PBL_DEPTH,        &
                       PSBL_DEPTH,PLMO,                              &
@@ -255,7 +250,7 @@ END MODULE MODI_TURB_VER
 !!           XCTV,XCHV   : cts for the T and moisture variances
 !!
 !!      Module MODD_PARAMETERS
-!!  
+!!
 !!           JPVEXT_TURB     : number of vertical external points
 !!           JPHEXT     : number of horizontal external points
 !!
@@ -310,6 +305,7 @@ END MODULE MODI_TURB_VER
 !!                                 reversed vertical levels
 !!                     10/2012 (J.Escobar) Bypass PGI bug , redefine some allocatable array inplace of automatic
 !!                     08/2014 (J.Escobar) Bypass PGI memory leak bug , replace IF statement with IF THEN ENDIF
+!!      Modifications: July,    2015  (Wim de Rooy) switch for HARATU (Racmo turbulence scheme)
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !! JL Redelsperger 03/2021 : add Ocean LES case
 !!--------------------------------------------------------------------------
@@ -317,10 +313,13 @@ END MODULE MODI_TURB_VER
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_DYN_n,          ONLY: LOCEAN
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_LES
@@ -339,11 +338,11 @@ USE MODI_TURB_VER_SV_FLUX
 USE MODI_TURB_VER_SV_CORR
 USE MODI_LES_MEAN_SUBGRID
 USE MODI_SBL_DEPTH
+USE MODI_SECOND_MNH
 !
 USE MODE_IO_FIELD_WRITE, only: IO_Field_write
 USE MODE_PRANDTL
 !
-USE MODI_SECOND_MNH
 !
 IMPLICIT NONE
 !
@@ -377,6 +376,9 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE    ! sinus of the angle
                                       ! between i and the slope vector
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
+! MFMOIST used in case of LHARATU
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST       ! moist mass flux dual scheme
+
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
                                                       ! Potential Temperature 
 !
@@ -405,6 +407,9 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PVSLOPEM     ! wind component along the
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
+! PLENGTHM PLENGTHH used in case of LHARATU
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLENGTHM     ! Turb. mixing length momentum
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLENGTHH     ! Turb. mixing length heat/moisture 
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
@@ -468,6 +473,7 @@ REAL, ALLOCATABLE, DIMENSION(:,:,:,:)  ::  &
        ZREDS1,   & ! 1D Redelsperger number R_sv
        ZRED2THS, & ! 3D Redelsperger number R*2_thsv
        ZRED2RS     ! 3D Redelsperger number R*2_rsv
+REAL, DIMENSION(SIZE(PLM,1),SIZE(PLM,2),SIZE(PLM,3))  ::  ZLM
 !
 LOGICAL :: GUSERV    ! flag to use water vapor
 INTEGER :: IKB,IKE   ! index value for the Beginning
@@ -475,6 +481,7 @@ INTEGER :: IKB,IKE   ! index value for the Beginning
 INTEGER :: JSV       ! loop counter on scalar variables
 REAL    :: ZTIME1
 REAL    :: ZTIME2
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
 TYPE(TFIELDDATA) :: TZFIELD
 !----------------------------------------------------------------------------
 ALLOCATE (      ZBETA(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3))    ,&
@@ -509,6 +516,8 @@ ALLOCATE ( &
 !*       1.   PRELIMINARIES
 !             -------------
 !
+IF (LHOOK) CALL DR_HOOK('TURB_VER',0,ZHOOK_HANDLE)
+!
 IKB=KKA+JPVEXT_TURB*KKL
 IKE=KKU-JPVEXT_TURB*KKL
 !
@@ -542,25 +551,26 @@ ZSQRT_TKE = SQRT(PTKEM)
 !
 ! gradients of mean quantities at previous time-step
 !
-ZDTH_DZ = GZ_M_W(KKA,KKU,KKL,PTHLM(:,:,:),PDZZ)
+ZDTH_DZ = GZ_M_W(PTHLM(:,:,:),PDZZ, KKA, KKU, KKL)
 ZDR_DZ  = 0.
-IF (KRR>0) THEN
-ZDR_DZ  = GZ_M_W(KKA,KKU,KKL,PRM(:,:,:,1),PDZZ)
-ENDIF
+IF (KRR>0) ZDR_DZ  = GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL)
 !
 !
 ! Denominator factor in 3rd order terms
 !
-ZD(:,:,:) = (1.+ZREDTH1+ZREDR1) * (1.+0.5*(ZREDTH1+ZREDR1))
+IF (.NOT. LHARAT) THEN
+  ZD(:,:,:) = (1.+ZREDTH1+ZREDR1) * (1.+0.5*(ZREDTH1+ZREDR1))
+ELSE
+  ZD(:,:,:) = 1.
+ENDIF
 !
 ! Phi3 and Psi3 Prandtl numbers
 !
 GUSERV = KRR/=0
 !
 ZPHI3 = PHI3(ZREDTH1,ZREDR1,ZRED2TH3,ZRED2R3,ZRED2THR3,HTURBDIM,GUSERV)
-IF(KRR/=0) THEN
+IF(KRR/=0) &
 ZPSI3 = PSI3(ZREDR1,ZREDTH1,ZRED2R3,ZRED2TH3,ZRED2THR3,HTURBDIM,GUSERV)
-ENDIF
 !
 ! Prandtl numbers for scalars
 !
@@ -591,6 +601,12 @@ END IF
 !*       4.   TURBULENT CORRELATIONS : <w Rc>, <THl THl>, <THl Rnp>, <Rnp Rnp>
 !             ----------------------------------------------------------------
 !
+
+IF (LHARAT) THEN
+  ZLM=PLENGTHH
+ELSE
+  ZLM=PLM
+ENDIF
 !
   CALL  TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI,               &
                         OTURB_FLX,HTURBDIM,HTOM,                      &
@@ -600,13 +616,14 @@ END IF
                         PRHODJ,PTHVREF,                               &
                         PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
                         PWM,PTHLM,PRM,PSVM,                           &
-                        PTKEM,PLM,PLEPS,                              &
+                        PTKEM,ZLM,PLEPS,                         &
                         PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE,    &
                         ZBETA, ZSQRT_TKE, ZDTH_DZ, ZDR_DZ, ZRED2TH3,  &
                         ZRED2R3, ZRED2THR3, ZBLL_O_E, ZETHETA,        &
                         ZEMOIST, ZREDTH1, ZREDR1, ZPHI3, ZPSI3, ZD,   &
-                        PFWTH,PFWR,PFTH2,PFR2,PFTHR,PBL_DEPTH,ZWTHV,  &
-                        PRTHLS,PRRS,ZTHLP,ZRP,PTP,PWTH,PWRC           )
+                        PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
+                        MFMOIST,PBL_DEPTH,ZWTHV,                      &
+                        PRTHLS,PRRS,ZTHLP,ZRP,PTP,PWTH,PWRC )
 !
   CALL  TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,               &
                         OTURB_FLX,HTURBDIM,HTOM,                      &
@@ -616,13 +633,13 @@ END IF
                         PRHODJ,PTHVREF,                               &
                         PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
                         PWM,PTHLM,PRM,PSVM,                           &
-                        PTKEM,PLM,PLEPS,                              &
+                        PTKEM,ZLM,PLEPS,                              &
                         PLOCPEXNM,PATHETA,PAMOIST,PSRCM,              &
                         ZBETA, ZSQRT_TKE, ZDTH_DZ, ZDR_DZ, ZRED2TH3,  &
                         ZRED2R3, ZRED2THR3, ZBLL_O_E, ZETHETA,        &
                         ZEMOIST, ZREDTH1, ZREDR1, ZPHI3, ZPSI3, ZD,   &
                         PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
-                        ZTHLP,ZRP,PSIGS                          )
+                        ZTHLP,ZRP,MFMOIST,PSIGS                  )
 !
 !----------------------------------------------------------------------------
 !
@@ -637,6 +654,9 @@ END IF
 !*       7.   DIAGNOSTIC COMPUTATION OF THE 1D <W W> VARIANCE
 !             -----------------------------------------------
 !
+!
+IF (LHARAT) ZLM=PLENGTHM
+!
 CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
                       OTURB_FLX,KRR,                                &
                       HTURBDIM,PIMPL,PEXPL,PTSTEP,                  &
@@ -646,7 +666,7 @@ CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
                       PRHODJ,                                       &
                       PCDUEFF,PTAU11M,PTAU12M,PTAU33M,              &
                       PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
-                      PTKEM,PLM,ZWU,ZWV,                            &
+                      PTKEM,ZLM,MFMOIST,ZWU,ZWV,                    &
                       PRUS,PRVS,PRWS,                               &
                       PDP,PTP                                       )
 !
@@ -656,6 +676,8 @@ CALL  TURB_VER_DYN_FLUX(KKA,KKU,KKL,                                &
 !*       8.   SOURCES OF PASSIVE SCALAR VARIABLES
 !             -----------------------------------
 !
+IF (LHARAT) ZLM=PLENGTHH
+!
 IF (SIZE(PSVM,4)>0)                                                 &
 CALL  TURB_VER_SV_FLUX(KKA,KKU,KKL,                                 &
                       OTURB_FLX,HTURBDIM,                           &
@@ -665,7 +687,7 @@ CALL  TURB_VER_SV_FLUX(KKA,KKU,KKL,                                 &
                       PRHODJ,PWM,                                   &
                       PSFSVM,PSFSVP,                                &
                       PSVM,                                         &
-                      PTKEM,PLM,ZPSI_SV,                            &
+                      PTKEM,ZLM,MFMOIST,ZPSI_SV,                    &
                       PRSVS,PWSV                                    )
 !
 !
@@ -675,7 +697,7 @@ CALL  TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,                   &
                       PTHLM,PRM,PTHVREF,                            &
                       PLOCPEXNM,PATHETA,PAMOIST,PSRCM,ZPHI3,ZPSI3,  &
                       PWM,PSVM,                                     &
-                      PTKEM,PLM,PLEPS,ZPSI_SV                       )
+                      PTKEM,ZLM,PLEPS,ZPSI_SV                       )
 !
 !
 !----------------------------------------------------------------------------
@@ -692,7 +714,7 @@ IF (SIZE(PSBL_DEPTH)>0) CALL SBL_DEPTH(IKB,IKE,PZZ,ZWU,ZWV,ZWTHV,PLMO,PSBL_DEPTH
 !             ------
 !
 !
-IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. .NOT. LHARAT) THEN
 !
 ! stores the Turbulent Prandtl number
 ! 
@@ -743,4 +765,5 @@ END IF
 !
 !
 !----------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('TURB_VER',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER                                                                                               
diff --git a/src/mesonh/turb/turb_ver_dyn_flux.f90 b/src/mesonh/turb/turb_ver_dyn_flux.f90
index 51bc4e7e1b868799b038e1e56b089bd2cb88ae22..b7c131ba1912c6ebc604dd502a04ab02bf33a49a 100644
--- a/src/mesonh/turb/turb_ver_dyn_flux.f90
+++ b/src/mesonh/turb/turb_ver_dyn_flux.f90
@@ -19,7 +19,7 @@ INTERFACE
                       PRHODJ,                                       &
                       PCDUEFF,PTAU11M,PTAU12M,PTAU33M,              &
                       PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
-                      PTKEM,PLM,PWU,PWV,                            &
+                      PTKEM,PLM,MFMOIST,PWU,PWV,                    &
                       PRUS,PRVS,PRWS,                               &
                       PDP,PTP                                       )
 !
@@ -277,6 +277,7 @@ END MODULE MODI_TURB_VER_DYN_FLUX
 !!                     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) LHARATU switch
 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !!      Q. Rodier      17/01/2019 : cleaning : remove cyclic conditions on DP and ZA
@@ -286,11 +287,14 @@ END MODULE MODI_TURB_VER_DYN_FLUX
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CONF
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_DYN_n,          ONLY: LOCEAN
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_LES
 USE MODD_NSV
@@ -305,7 +309,8 @@ USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SECOND_MNH
-USE MODI_SHUMAN
+USE MODI_SHUMAN , ONLY: MZM, MZF, MXM, MXF, MYM, MYF,&
+                      & DZM, DXF, DXM, DYF, DYM
 USE MODI_TRIDIAG 
 USE MODI_TRIDIAG_WIND 
 USE MODI_LES_MEAN_SUBGRID
@@ -342,6 +347,8 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE    ! sinus of the angle
                                       ! between i and the slope vector
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
+
 !
 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PCDUEFF     ! Cd * || u || at time t
 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTAU11M      ! <uu> in the axes linked 
@@ -767,7 +774,7 @@ IF (LOCEAN) THEN
   ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE) 
 END IF
 !
-IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
   ! stores the V wind component vertical flux
   TZFIELD%CMNHNAME   = 'VW_VFLX'
   TZFIELD%CSTDNAME   = ''
@@ -900,7 +907,7 @@ END IF
 !*       7.   DIAGNOSTIC COMPUTATION OF THE 1D <W W> VARIANCE
 !             -----------------------------------------------
 !
-IF ( OTURB_FLX .AND. tpfile%lopened .AND. HTURBDIM == '1DIM') THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. HTURBDIM == '1DIM') THEN
   ZFLXZ(:,:,:)= (2./3.) * PTKEM(:,:,:)                     &
      -XCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ)
   ! to be tested &
@@ -921,4 +928,5 @@ END IF
 !
 !----------------------------------------------------------------------------
 !
+IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_DYN_FLUX
diff --git a/src/mesonh/turb/turb_ver_sv_corr.f90 b/src/mesonh/turb/turb_ver_sv_corr.f90
index b62268e7e82a28d876844fb7abbbcaa02f3e87d0..6676e227c4617cfe1ed9752796bd296d7b15abe4 100644
--- a/src/mesonh/turb/turb_ver_sv_corr.f90
+++ b/src/mesonh/turb/turb_ver_sv_corr.f90
@@ -97,6 +97,9 @@ END MODULE MODI_TURB_VER_SV_CORR
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_PARAMETERS
@@ -110,7 +113,7 @@ USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
-USE MODI_SHUMAN 
+USE MODI_SHUMAN , ONLY : MZF
 USE MODI_EMOIST
 USE MODI_ETHETA
 USE MODI_LES_MEAN_SUBGRID
@@ -168,6 +171,8 @@ REAL :: ZCTSVD = 2.4  ! constant for temperature - scalar covariance dissipation
 REAL :: ZCQSVD = 2.4  ! constant for humidity - scalar covariance dissipation
 !----------------------------------------------------------------------------
 !
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',0,ZHOOK_HANDLE)
 CALL SECOND_MNH(ZTIME1)
 !
 IF(LBLOWSNOW) THEN
@@ -220,4 +225,5 @@ CALL SECOND_MNH(ZTIME2)
 XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
 !----------------------------------------------------------------------------
 !
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_SV_CORR
diff --git a/src/mesonh/turb/turb_ver_sv_flux.f90 b/src/mesonh/turb/turb_ver_sv_flux.f90
index 23d8bee0342d4fc3a6f28f73733d8f35c27308bc..aa7010e432a73086a2e918a15d98bab9602c913a 100644
--- a/src/mesonh/turb/turb_ver_sv_flux.f90
+++ b/src/mesonh/turb/turb_ver_sv_flux.f90
@@ -18,7 +18,7 @@ INTERFACE
                       PRHODJ,PWM,                                   &
                       PSFSVM,PSFSVP,                                &
                       PSVM,                                         &
-                      PTKEM,PLM,PPSI_SV,                            &
+                      PTKEM,PLM,MFMOIST,PPSI_SV,                    &
                       PRSVS,PWSV                                    )
 !
 USE MODD_IO, ONLY: TFILEDATA
@@ -259,6 +259,7 @@ END MODULE MODI_TURB_VER_SV_FLUX
 !!                                              change of YCOMMENT
 !!                     Feb 2012(Y. Seity) add possibility to run with reversed 
 !!                                              vertical levels
+!!      Modifications: July 2015 (Wim de Rooy) LHARAT switch
 !!                     Feb 2017(M. Leriche) add initialisation of ZSOURCE
 !!                                   to avoid unknwon values outside physical domain
 !!                                   and avoid negative values in sv tendencies
@@ -268,22 +269,25 @@ END MODULE MODI_TURB_VER_SV_FLUX
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_LES
 USE MODD_CONF
 USE MODD_NSV,            ONLY: XSVMIN, NSV_LGBEG, NSV_LGEND
 USE MODD_BLOWSNOW
-USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 !
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
-USE MODI_SHUMAN 
+USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
 USE MODI_TRIDIAG 
 USE MODI_TRIDIAG_WIND 
 USE MODI_EMOIST
@@ -314,6 +318,8 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
                                                       ! normal to the ground surface
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST       ! moist mf dual scheme
+
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSFSVM       ! t - deltat 
 !
@@ -365,6 +371,9 @@ TYPE(TFIELDDATA)  :: TZFIELD
 !*       1.   PRELIMINARIES
 !             -------------
 !
+
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',0,ZHOOK_HANDLE)
 IKB=KKA+JPVEXT_TURB*KKL
 IKE=KKU-JPVEXT_TURB*KKL
 IKT=SIZE(PSVM,3)
@@ -424,7 +433,7 @@ DO JSV=1,ISV
 ! PRSVS(:,:,:,JSV)= MAX((PRSVS(:,:,:,JSV)+    &
 !                   PRHODJ(:,:,:)*(ZRES(:,:,:)-PSVM(:,:,:,JSV))/PTSTEP),XSVMIN(JSV))
 !
-  IF ( (OTURB_FLX .AND. tpfile%lopened) .OR. LLES_CALL ) THEN
+  IF ( (OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL ) THEN
     ! Diagnostic of the cartesian vertical flux
     !
     ZFLXZ(:,:,:) = -ZCSV * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM)) / PDZZ * &
@@ -451,7 +460,7 @@ DO JSV=1,ISV
     PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-KKL,JSV)
  END IF
   !
-  IF (OTURB_FLX .AND. tpfile%lopened) THEN
+  IF (OTURB_FLX .AND. TPFILE%LOPENED) THEN
     ! stores the JSVth vertical flux
     WRITE(TZFIELD%CMNHNAME,'("WSV_FLX_",I3.3)') JSV
     TZFIELD%CSTDNAME   = ''
@@ -487,4 +496,5 @@ END DO   ! end of scalar loop
 !
 !----------------------------------------------------------------------------
 !
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_SV_FLUX
diff --git a/src/mesonh/turb/turb_ver_thermo_corr.f90 b/src/mesonh/turb/turb_ver_thermo_corr.f90
index bdd074e5c52af78809d84b8fa5077d56e12a76d5..6edf5c72431dddc3b9aeb587a09d7a8bd1a5f45b 100644
--- a/src/mesonh/turb/turb_ver_thermo_corr.f90
+++ b/src/mesonh/turb/turb_ver_thermo_corr.f90
@@ -23,106 +23,7 @@ INTERFACE
                       PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
                       PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
                       PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
-                      PTHLP,PRP,PSIGS                          )
-!
-USE MODD_IO, ONLY: TFILEDATA
-!
-INTEGER,                INTENT(IN)   :: KKA           !near ground array index  
-INTEGER,                INTENT(IN)   :: KKU           !uppest atmosphere array index
-INTEGER,                INTENT(IN)   :: KKL           !vert. levels type 1=MNH -1=AR 
-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.
-LOGICAL,                INTENT(IN)   ::  OTURB_FLX    ! switch to write the
-                                 ! turbulent fluxes in the syncronous FM-file
-CHARACTER(len=4),       INTENT(IN)   ::  HTURBDIM     ! dimensionality of the
-                                                      ! turbulence scheme
-CHARACTER(len=4),       INTENT(IN)   ::  HTOM         ! type of Third Order Moment
-REAL,                   INTENT(IN)   ::  PIMPL, PEXPL ! Coef. for temporal disc.
-TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDZZ, PDXX, PDYY, PDZX, PDZY
-                                                      ! Metric coefficients
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
-                                                      ! normal to the ground surface
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
-                                                      ! Potential Temperature 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHM,PSFRM ! surface fluxes at time
-!                                                     ! t - deltat 
-!
-REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTHP,PSFRP ! surface fluxes at time
-!                                                     ! t + deltat 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PWM 
-! Vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHLM 
-! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios 
-                                                      ! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! Mixing ratios 
-!
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PATHETA      ! coefficients between 
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PAMOIST      ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSRCM        ! normalized 
-                  ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBETA        ! buoyancy coefficient
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PSQRT_TKE    ! sqrt(e)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDTH_DZ      ! d(th)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PDR_DZ       ! d(rt)/dz
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2TH3     ! 3D Redeslperger number R*2_th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2R3      ! 3D Redeslperger number R*2_r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRED2THR3    ! 3D Redeslperger number R*2_thr
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PBLL_O_E     ! beta * Lk * Leps / tke
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PETHETA      ! Coefficient for theta in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PEMOIST      ! Coefficient for r in theta_v computation
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDTH1      ! 1D Redelsperger number for Th
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PREDR1       ! 1D Redelsperger number for r
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPHI3        ! Prandtl number for temperature
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PPSI3        ! Prandtl number for vapor
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PD           ! Denominator in Prandtl numbers
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWTH        ! d(w'2th' )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFWR         ! d(w'2r'  )/dz (at flux point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTH2        ! d(w'th'2 )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFR2         ! d(w'r'2  )/dz (at mass point)
-REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PFTHR        ! d(w'th'r')/dz (at mass point)
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHLP      ! guess of thl at t+ deltat
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRP        ! guess of r at t+ deltat
-!
-REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PSIGS     ! Vert. part of Sigma_s at t
-!
-!
-!
-END SUBROUTINE TURB_VER_THERMO_CORR
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_VER_THERMO_CORR
-!
-!
-!     ###############################################################
-      SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR, KRRL, KRRI,  &
-                      OTURB_FLX,HTURBDIM,HTOM,                      &
-                      PIMPL,PEXPL,                                  &
-                      TPFILE,                                       &
-                      PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,           &
-                      PRHODJ,PTHVREF,                               &
-                      PSFTHM,PSFRM,PSFTHP,PSFRP,                    &
-                      PWM,PTHLM,PRM,PSVM,                           &
-                      PTKEM,PLM,PLEPS,                              &
-                      PLOCPEXNM,PATHETA,PAMOIST,PSRCM,              &
-                      PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3,  &
-                      PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
-                      PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
-                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,                  &
-                      PTHLP,PRP,PSIGS                          )
+                      PTHLP,PRP,MFMOIST,PSIGS                       )
 !     ###############################################################
 !
 !
@@ -304,15 +205,18 @@ END MODULE MODI_TURB_VER_THERMO_CORR
 !!                                              change of YCOMMENT
 !!                     2012-02 (Y. Seity) add possibility to run with reversed 
 !!                                              vertical levels
+!!      Modifications  July 2015 (Wim de Rooy) LHARAT switch
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 USE MODD_CST
 USE MODD_CTURB
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS
 USE MODD_CONF
@@ -322,7 +226,7 @@ USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
-USE MODI_SHUMAN 
+USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
 USE MODI_TRIDIAG 
 USE MODI_LES_MEAN_SUBGRID
 USE MODI_PRANDTL
@@ -359,6 +263,8 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
                                                       ! normal to the ground surface
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
+
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
                                                       ! Potential Temperature 
 !
@@ -377,6 +283,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios
 REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! Mixing ratios 
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
+! In case LHARATU=TRUE, PLM already includes all stability corrections
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
@@ -567,7 +474,7 @@ END IF
   !
   !
   ! stores <THl THl>  
-  IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
     TZFIELD%CMNHNAME   = 'THL_VVAR'
     TZFIELD%CSTDNAME   = ''
     TZFIELD%CLONGNAME  = 'THL_VVAR'
@@ -694,7 +601,7 @@ END IF
                      2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLXZ(:,:,:)
     END IF
     ! stores <THl Rnp>   
-    IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+    IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
       TZFIELD%CMNHNAME   = 'THLRCONS_VCOR'
       TZFIELD%CSTDNAME   = ''
       TZFIELD%CLONGNAME  = 'THLRCONS_VCOR'
@@ -801,7 +708,7 @@ END IF
       PSIGS(:,:,:) = PSIGS(:,:,:) + PAMOIST(:,:,:) **2 * ZFLXZ(:,:,:)
     END IF
     ! stores <Rnp Rnp>    
-    IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+    IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
       TZFIELD%CMNHNAME   = 'RTOT_VVAR'
       TZFIELD%CSTDNAME   = ''
       TZFIELD%CLONGNAME  = 'RTOT_VVAR'
@@ -845,4 +752,5 @@ END IF
 !
   DEALLOCATE(ZCOEFF)
 !----------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_CORR',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_THERMO_CORR
diff --git a/src/mesonh/turb/turb_ver_thermo_flux.f90 b/src/mesonh/turb/turb_ver_thermo_flux.f90
index cf539984e7751f862a6251a3e80b048274740073..698bcfa762a20dd136103f3dc4793da6cd41a4c6 100644
--- a/src/mesonh/turb/turb_ver_thermo_flux.f90
+++ b/src/mesonh/turb/turb_ver_thermo_flux.f90
@@ -23,7 +23,7 @@ INTERFACE
                       PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3,  &
                       PRED2R3, PRED2THR3, PBLL_O_E, PETHETA,        &
                       PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD,   &
-                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,PBL_DEPTH,        &
+                      PFWTH,PFWR,PFTH2,PFR2,PFTHR,MFMOIST,PBL_DEPTH,&
                       PWTHV,PRTHLS,PRRS,PTHLP,PRP,PTP,PWTH,PWRC     )
 !
 USE MODD_IO, ONLY: TFILEDATA
@@ -322,6 +322,7 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
 !!                                              change of YCOMMENT
 !!                     2012-02 (Y. Seity) add possibility to run with reversed
 !!                                             vertical levels
+!!      Modifications  July 2015 (Wim de Rooy) LHARAT switch
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !!                     2021 (D. Ricard) last version of HGRAD turbulence scheme
 !!                                 Leronard terms instead of Reynolds terms
@@ -336,9 +337,12 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 USE MODD_CST
 USE MODD_CTURB
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD,          ONLY: TFIELDDATA, TYPEREAL
 USE MODD_GRID_n,         ONLY: XZS, XXHAT, XYHAT
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_METRICS_n,      ONLY: XDXX, XDYY, XDZX, XDZY, XDZZ
@@ -367,7 +371,7 @@ USE MODI_PRANDTL
 USE MODI_TRIDIAG_THERMO
 USE MODI_TM06_H
 !
-USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_PRANDTL
 !
 USE MODI_SECOND_MNH
@@ -402,6 +406,7 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PDIRCOSZW    ! Director Cosinus of the
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PZZ          ! altitudes
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PRHODJ       ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  MFMOIST      ! moist mass flux dual scheme
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHVREF      ! ref. state Virtual 
                                                       ! Potential Temperature 
 !
@@ -420,6 +425,8 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM          ! Mixing ratios
 REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PSVM         ! Mixing ratios 
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTKEM        ! TKE at time t
+!
+! In case LHARAT=TRUE, PLM already includes all stability corrections
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLM          ! Turb. mixing length   
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLEPS        ! dissipative length   
 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PLOCPEXNM    ! Lv(T)/Cp/Exnref at time t-1
@@ -751,7 +758,7 @@ ELSE
   PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
 END IF
 !
-IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
   ! stores the conservative potential temperature vertical flux
   TZFIELD%CMNHNAME   = 'THW_FLX'
   TZFIELD%CSTDNAME   = ''
@@ -987,7 +994,7 @@ IF (KRR /= 0) THEN
   PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
   !
   !
-  IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
     ! stores the conservative mixing ratio vertical flux
     TZFIELD%CMNHNAME   = 'RCONSW_FLX'
     TZFIELD%CSTDNAME   = ''
@@ -1063,7 +1070,7 @@ END IF
 !
 !*       4.1  <w Rc>    
 !
-IF ( ((OTURB_FLX .AND. tpfile%lopened) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
+IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
   !  
   ! recover the Conservative potential temperature flux : 
   ZA(:,:,:)   = DZM(PIMPL * PTHLP + PEXPL * PTHLM) / PDZZ *       &
@@ -1077,7 +1084,7 @@ IF ( ((OTURB_FLX .AND. tpfile%lopened) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
   ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB) 
   !                 
   ! store the liquid water mixing ratio vertical flux
-  IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
+  IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
     TZFIELD%CMNHNAME   = 'RCW_FLX'
     TZFIELD%CSTDNAME   = ''
     TZFIELD%CLONGNAME  = 'RCW_FLX'
@@ -1106,4 +1113,5 @@ IF (LOCEAN.AND.LDEEPOC) THEN
 END IF
 !
 !----------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',1,ZHOOK_HANDLE)
 END SUBROUTINE TURB_VER_THERMO_FLUX