diff --git a/docs/Interfaces b/docs/Interfaces
new file mode 100644
index 0000000000000000000000000000000000000000..a233c0e9d195229feaf8887c232326d1b096a995
--- /dev/null
+++ b/docs/Interfaces
@@ -0,0 +1,25 @@
+The PHYEX parameterizations can be called from different models (eg. Meso-NH, AROME) and from
+a driver (which will be included in this repository).
+Moreover, PHYEX parameterizations call externals subroutines, which are dependencies.
+
+This document aims at listing the interfaces to call the PHYEX parameterizations and the
+interfaces of the external modules called by PHYEX.
+
+PHYEX interfaces:
+- lima_adjust
+- ice_adjust
+- shallow_mf
+- turb
+- lima, lima_warm, lima_cold et lima_mixed
+- rain_ice, rain_ice_old
+- ini_...
+
+Dependencies:
+- budget
+- mode_msg, modd_io
+- modd_precision
+- modd_cst
+- yomhook, parkind1
+
+
+
diff --git a/docs/TODO b/docs/TODO
new file mode 100644
index 0000000000000000000000000000000000000000..eb46ff8bb1c42921906761f37106db3184de1c94
--- /dev/null
+++ b/docs/TODO
@@ -0,0 +1,16 @@
+LOCEAN:
+  La clé LOCEAN est dans un module spécifique à Méso-NH (MODD_DYNn).
+  Une solution serait de créer un module propre à PHYEX qui contiendrait des clés de contrôle de haut niveau
+  pour la physique (y en a-t-il d'autres?). Ce module serait initialisé dans Méso-NH à partir de la clé
+  actuelle qui est sans doute utilisée ailleurs dans le code de Méso-NH
+
+Dependencies:
+  - définir les interfaces propres
+  - créer des codes pour le driver
+  - liste dans document Interfaces
+  - pour AROME placés, en attendant, dans phyex/externals
+
+Clé de compilation REPRO48 ajoutée pour permettre de reproduire le cycle 48, elle:
+- contourne des corrections de bug
+- modifie l'organisation de calculs
+Cette clé devra être supprimée
diff --git a/src/arome/micro/budget.F90 b/src/arome/externals/budget.F90
similarity index 100%
rename from src/arome/micro/budget.F90
rename to src/arome/externals/budget.F90
diff --git a/src/arome/externals/modd_io.F90 b/src/arome/externals/modd_io.F90
new file mode 100644
index 0000000000000000000000000000000000000000..c111e469fb3d2014b7021eb698d564a94f0cf413
--- /dev/null
+++ b/src/arome/externals/modd_io.F90
@@ -0,0 +1,6 @@
+MODULE MODD_IO
+IMPLICIT NONE
+
+INTEGER, PARAMETER :: NVERB_NO=0, NVERB_FATAL=1, NVERB_ERROR=2, NVERB_WARNING=3, NVERB_INFO=4, NVERB_DEBUG=5
+INTEGER, SAVE :: N_ABORT_LEVEL = NVERB_ERROR
+ENDMODULE MODD_IO
diff --git a/src/arome/externals/modd_precision.F90 b/src/arome/externals/modd_precision.F90
new file mode 100644
index 0000000000000000000000000000000000000000..09509322d90cfeb2ace272ac779a2a7fdc582c09
--- /dev/null
+++ b/src/arome/externals/modd_precision.F90
@@ -0,0 +1,8 @@
+MODULE MODD_PRECISION
+USE PARKIND1
+IMPLICIT NONE
+SAVE
+
+INTEGER, PARAMETER :: MNHREAL = JPRB
+INTEGER, PARAMETER :: MNHREAL64 = JPRD
+ENDMODULE MODD_PRECISION
diff --git a/src/arome/externals/mode_msg.F90 b/src/arome/externals/mode_msg.F90
new file mode 100644
index 0000000000000000000000000000000000000000..d6276382965b1b5b1bc27ea2d5314cc8f7547210
--- /dev/null
+++ b/src/arome/externals/mode_msg.F90
@@ -0,0 +1,115 @@
+! Author(s)
+!   S. Riette (18 Nov 2021), adapted from the Meso-NH version
+! Modifications:
+!-----------------------------------------------------------------
+MODULE MODE_MSG
+
+USE MODD_IO, ONLY: NVERB_FATAL, NVERB_ERROR, NVERB_WARNING, &
+                  &NVERB_INFO, NVERB_DEBUG, N_ABORT_LEVEL
+
+IMPLICIT NONE
+
+INTEGER, PARAMETER :: NMSGLGTMAX   = 100 ! Maximum length for a message
+INTEGER, PARAMETER :: NMSGLLINEMAX = 10  ! Maximum number of lines for a message
+CHARACTER(LEN=NMSGLGTMAX), DIMENSION(NMSGLLINEMAX) :: CMNHMSG=''
+
+#include "abor1.intfb.h"
+
+INTERFACE PRINT_MSG
+  MODULE PROCEDURE PRINT_MSG_1LINE, PRINT_MSG_MULTI_CMNHMSG, PRINT_MSG_MULTI
+ENDINTERFACE PRINT_MSG
+
+CONTAINS
+
+SUBROUTINE PRINT_MSG_1LINE(KVERB, HDOMAIN, HSUBR, HMSG)
+  INTEGER,          INTENT(IN) :: KVERB   !Verbosity level
+  CHARACTER(LEN=*), INTENT(IN) :: HDOMAIN !Domain/category of message
+  CHARACTER(LEN=*), INTENT(IN) :: HSUBR   !Subroutine/function name
+  CHARACTER(LEN=*), INTENT(IN) :: HMSG    !Message
+
+  CALL PRINT_MSG_MULTI(KVERB, HDOMAIN, HSUBR, [HMSG])
+
+ENDSUBROUTINE PRINT_MSG_1LINE
+
+SUBROUTINE PRINT_MSG_MULTI_CMNHMSG(KVERB, HDOMAIN, HSUBR)
+  INTEGER,          INTENT(IN) :: KVERB   !Verbosity level
+  CHARACTER(LEN=*), INTENT(IN) :: HDOMAIN !Domain/category of message
+  CHARACTER(LEN=*), INTENT(IN) :: HSUBR   !Subroutine/function name
+
+  INTEGER :: ILINES
+
+  !Find the last non empty line
+  ILINES=SIZE(CMNHMSG)
+  DO WHILE (LEN_TRIM(CMNHMSG(ILINES))==0)
+    ILINES=ILINES - 1
+  ENDDO
+
+  CALL PRINT_MSG_MULTI(KVERB, HDOMAIN, HSUBR, CMNHMSG(1:ILINES))
+
+  !Empty the message buffer
+  !This is necessary especially if the next call contain a shorter message
+  CMNHMSG(1:ILINES)=''
+
+ENDSUBROUTINE PRINT_MSG_MULTI_CMNHMSG
+
+SUBROUTINE PRINT_MSG_MULTI(KVERB, HDOMAIN, HSUBR, HMSG)
+!
+USE YOMLUN, ONLY : NULOUT
+!
+!
+INTEGER,                        INTENT(IN) :: KVERB   !Verbosity level
+CHARACTER(LEN=*),               INTENT(IN) :: HDOMAIN !Domain/category of message
+CHARACTER(LEN=*),               INTENT(IN) :: HSUBR   !Subroutine/function name
+CHARACTER(LEN=*), dimension(:), INTENT(IN) :: HMSG    !Message
+!
+CHARACTER(LEN=2)  :: YSZ
+CHARACTER(LEN=9)  :: YPRE
+CHARACTER(LEN=30) :: YSUBR
+CHARACTER(LEN=:), ALLOCATABLE :: YFORMAT
+INTEGER :: JI
+INTEGER :: ILINES
+!
+ILINES=SIZE(HMSG)
+
+SELECT CASE(KVERB)
+  CASE(NVERB_FATAL)
+    YPRE='FATAL:   '
+  CASE(NVERB_ERROR)
+    YPRE='ERROR:   '
+  CASE(NVERB_WARNING)
+    YPRE='WARNING: '
+  CASE(NVERB_INFO)
+    YPRE='INFO:    '
+  CASE(NVERB_DEBUG)
+    YPRE='DEBUG:   '
+  CASE DEFAULT
+    WRITE(UNIT=NULOUT, FMT=*) 'ERROR: PRINT_MSG: wrong verbosity level'
+END SELECT
+!
+YSUBR=TRIM(HSUBR)//':'
+
+IF (ILINES==1) THEN
+  WRITE(UNIT=NULOUT, FMT="(A9,A30,A)") YPRE, YSUBR, TRIM(HMSG(1))
+ELSE
+ IF (ILINES<10) THEN
+    YSZ = 'I1'
+  ELSEIF (ILINES<100) THEN
+    YSZ = 'I2'
+  ELSEIF (ILINES<1000) THEN
+    YSZ = 'I3'
+  ELSE
+    YSZ = 'I4'
+  ENDIF
+  YFORMAT='(A9,A30,' // YSZ // ',''/'',' // YSZ // ','': '',A)'
+  DO JI=1, ILINES
+    WRITE(UNIT=NULOUT, FMT=YFORMAT) YPRE, YSUBR, JI, ILINES, TRIM(HMSG(JI))
+  ENDDO
+ENDIF
+!
+IF (KVERB<=N_ABORT_LEVEL) THEN
+  CALL ABOR1(TRIM(HMSG(ILINES))) !Last line repeated
+END IF
+!
+ENDSUBROUTINE PRINT_MSG_MULTI
+
+ENDMODULE MODE_MSG
diff --git a/src/arome/gmkpack_ignored_files b/src/arome/gmkpack_ignored_files
new file mode 100644
index 0000000000000000000000000000000000000000..e45272470ccf31f649b9aac1c2ae0c19551d472c
--- /dev/null
+++ b/src/arome/gmkpack_ignored_files
@@ -0,0 +1 @@
+phyex/micro/budget.F90
diff --git a/src/arome/micro/ini_cst.F90 b/src/arome/micro/ini_cst.F90
index 82d60bc694f7c25c8cbf2c2615f6ea7248a5ef72..f0f6b11870c2ddd81f19531c01ee76719870a13f 100644
--- a/src/arome/micro/ini_cst.F90
+++ b/src/arome/micro/ini_cst.F90
@@ -94,6 +94,13 @@ XG      = 9.80665
 !*	 4.     REFERENCE PRESSURE
 !	        -------------------
 !
+! Ocean model cst same as in 1D/CMO SURFEX
+! values used in ini_cst to overwrite XP00 and XTH00
+XRH00OCEAN =1024.
+XTH00OCEAN = 286.65
+XSA00OCEAN= 32.6
+XP00OCEAN = 201.E5
+!Atmospheric model
 XP00 = 1.E5
 XTH00 = 300.
 !-------------------------------------------------------------------------------
@@ -133,7 +140,11 @@ XALPW  = LOG(XESTT) + (XBETAW /XTT) + (XGAMW *LOG(XTT))
 XGAMI  = (XCI - XCPV) / XRV
 XBETAI = (XLSTT/XRV) + (XGAMI * XTT)
 XALPI  = LOG(XESTT) + (XBETAI /XTT) + (XGAMI *LOG(XTT))
-!-------------------------------------------------------------------------------
+! Values identical to ones used in CMO1D in SURFEX /could be modified
+! Coefficient of thermal expansion of water (K-1)
+XALPHAOC = 1.9E-4
+! Coeff of Haline contraction coeff (S-1)
+XBETAOC= 7.7475E-4
 !
 !*	 7.     PRECOMPUTED CONSTANTS
 !	        ---------------------
diff --git a/src/arome/micro/modd_cst.F90 b/src/arome/micro/modd_cst.F90
index 806be1b8441c277dcd29822149f4b2e7deaa699a..1f5d39b521df169340f648dcd59e73afbe32b217 100644
--- a/src/arome/micro/modd_cst.F90
+++ b/src/arome/micro/modd_cst.F90
@@ -32,6 +32,7 @@
 !!      C. Mari     31/10/00  add NDAYSEC
 !!      V. Masson   01/03/03  add conductivity of ice
 !!      R. El Khatib 04/08/14 add pre-computed quantities
+!!      J.L. Redelsperger 03/2021  add constants for ocean penetrating solar
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -53,12 +54,14 @@ REAL,SAVE :: XRADIUS,XOMEGA     ! Earth radius, earth rotation
 REAL,SAVE :: XG                 ! Gravity constant
 !
 REAL,SAVE :: XP00               ! Reference pressure
+REAL,SAVE :: XP00OCEAN          ! Reference pressure for ocean model
+REAL,SAVE :: XRH00OCEAN         ! Reference density for ocean model
 !
 REAL,SAVE :: XSTEFAN,XI0        ! Stefan-Boltzman constant, solar constant
 !
 REAL,SAVE :: XMD,XMV            ! Molar mass of dry air and molar mass of vapor
 REAL,SAVE :: XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
-REAL,SAVE :: XEPSILO            ! XMV/XMD 
+REAL,SAVE :: XEPSILO            ! XMV/XMD
 REAL,SAVE :: XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
 REAL,SAVE :: XRHOLW             ! Volumic mass of liquid water
 REAL,SAVE :: XCL,XCI            ! Cl (liquid), Ci (ice)
@@ -73,8 +76,19 @@ REAL,SAVE :: XALPW,XBETAW,XGAMW ! Constants for saturation vapor
 REAL,SAVE :: XALPI,XBETAI,XGAMI ! Constants for saturation vapor
                                 !  pressure  function over solid ice
 REAL,SAVE :: XCONDI             ! thermal conductivity of ice (W m-1 K-1)
-REAL, SAVE        :: XTH00      ! reference value  for the potential
-                                ! temperature
+REAL,SAVE :: XALPHAOC           ! thermal expansion coefficient for ocean (K-1)
+REAL,SAVE :: XBETAOC             ! Haline contraction coeff for ocean (S-1)
+REAL,SAVE :: XTH00              ! reference value  for the potential temperature
+REAL,SAVE :: XTH00OCEAN         ! Ref value for pot temp in ocean model
+REAL,SAVE :: XSA00OCEAN         ! Ref value for SAlinity in ocean model
+REAL,SAVE :: XROC=0.69! 3 coeffs for SW penetration in  Ocean (Hoecker et al)
+REAL,SAVE :: XD1=1.1
+REAL,SAVE :: XD2=23.
+! Values used in SURFEX CMO
+!REAL,SAVE :: XROC=0.58
+!REAL,SAVE :: XD1=0.35
+!REAL,SAVE :: XD2=23.
+
 REAL,SAVE :: XRHOLI             ! Volumic mass of liquid water
 !
 INTEGER, SAVE :: NDAYSEC        ! Number of seconds in a day
diff --git a/src/arome/turb/compute_updraft.F90 b/src/arome/turb/compute_updraft.F90
index 3480e0d50915cbd66658b2db79520cca366f670b..8ba4d1c567308de0dde16fbc00d96e017c4c4a4a 100644
--- a/src/arome/turb/compute_updraft.F90
+++ b/src/arome/turb/compute_updraft.F90
@@ -49,6 +49,7 @@
 !!     S. Riette Jan 2012: support for both order of vertical levels
 !!     V.Masson, C.Lac : 02/2011 : SV_UP initialized by a non-zero value
 !!     S. Riette Apr 2013: improvement of continuity at the condensation level
+!!     Q.Rodier  01/2019 : support RM17 mixing length
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -56,6 +57,7 @@
 !
 USE MODD_CST
 USE MODD_CMFSHALL
+USE MODD_TURB_n, ONLY : CTURBLEN
 
 USE MODI_COMPUTE_ENTR_DETR
 USE MODI_TH_R_FROM_THL_RT_1D
@@ -170,6 +172,7 @@ REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP,&
 REAL  :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
 
 REAL  :: ZTMAX,ZRMAX  ! control value
+REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT',0,ZHOOK_HANDLE)
 
@@ -291,8 +294,19 @@ IF (OENTR_DETR) THEN
   ! compute L_up
   GLMIX=.TRUE.
   ZTKEM_F(:,KKB)=0.
-
-  CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.TRUE.,ZLUP)
+  !
+  IF(CTURBLEN=='RM17') THEN
+    ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
+    ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
+    ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+    PRINT*, 'phasage bete sans controle'
+    CALL ABORT
+    STOP
+  ELSE
+    ZSHEAR = 0. !no shear in bl89 mixing length
+  END IF
+ !
+  CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
   ZLUP(:)=MAX(ZLUP(:),1.E-10)
 
   ! Compute Buoyancy flux at the ground
diff --git a/src/arome/turb/compute_updraft_rhcj10.F90 b/src/arome/turb/compute_updraft_rhcj10.F90
index e28b53ad9a2fa9b4083fed7da802a1315c729962..14f0b6ae68f2f53203e9096e4cdb74f6bae6fc48 100644
--- a/src/arome/turb/compute_updraft_rhcj10.F90
+++ b/src/arome/turb/compute_updraft_rhcj10.F90
@@ -44,6 +44,7 @@
 !!     Y. Bouteloup (2012)
 !!     R. Honert Janv 2013 ==> corection of some coding bugs
 !!     R. El Khatib 15-Oct-2014 Optimization
+!!     Q.Rodier  01/2019 : support RM17 mixing length
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -51,7 +52,7 @@
 !
 USE MODD_CST
 USE MODD_CMFSHALL
-
+USE MODD_TURB_n, ONLY : CTURBLEN
 USE MODI_TH_R_FROM_THL_RT_1D
 USE MODI_SHUMAN_MF
 
@@ -177,6 +178,7 @@ REAL  :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
 
 REAL  :: ZTMAX,ZRMAX, ZEPS  ! control value
 
+REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT_RHCJ10',0,ZHOOK_HANDLE)
 
@@ -322,9 +324,20 @@ ENDDO
 ! compute L_up
 GLMIX=.TRUE.
 ZTKEM_F(:,KKB)=0.
-
-
-CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM_F,KKB,GLMIX,.FALSE.,ZLUP)
+!
+IF(CTURBLEN=='RM17') THEN
+  ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
+  ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
+  ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
+  PRINT*, 'phasage bete sans controle'
+  CALL ABORT
+  STOP
+ELSE
+  ZSHEAR = 0. !no shear in bl89 mixing length
+END IF
+!
+CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
+                       ZTHVM_F,KKB,GLMIX,.FALSE.,ZSHEAR,ZLUP)
 ZLUP(:)=MAX(ZLUP(:),1.E-10)
 
 DO JI=1,IIJU
diff --git a/src/arome/turb/ini_cturb.F90 b/src/arome/turb/ini_cturb.F90
index dec13f81fac495967bd53b6a0b88a804876cc061..06bbc545492660b517fea9f828a6abd13ee66937 100644
--- a/src/arome/turb/ini_cturb.F90
+++ b/src/arome/turb/ini_cturb.F90
@@ -144,6 +144,10 @@ XALPSBL = 4.63
 !       Stull 1988                  = 4.75
 !
 !
+!         1.11  Value related to the shear term in mixing length computation
+!
+XRM17 = 0.5  ! Rodier et al 2017
+!
 !
 !         2. Derivated constants
 !            -------------------
diff --git a/src/arome/turb/modd_cturb.F90 b/src/arome/turb/modd_cturb.F90
index 49ecfb6777011f51d98bdda53bec3ad5926ca7a2..9b18f803fd4181c9136220c02db637666cc20d0a 100644
--- a/src/arome/turb/modd_cturb.F90
+++ b/src/arome/turb/modd_cturb.F90
@@ -58,6 +58,7 @@ REAL,SAVE :: XCDD         ! ct. for the destruction term in the dissipation eq.
 REAL,SAVE :: XCDT         ! ct. for the transport term in the dissipation eq.
 !
 REAL,SAVE :: XTKEMIN      ! mimimum value for the TKE
+REAL,SAVE :: XRM17        ! Rodier et al 2017 constant in shear term for mixing length
 !
 REAL,SAVE :: XLINI        ! initial value for BL mixing length
 REAL,SAVE :: XLINF        ! to prevent division by zero in the BL algorithm
diff --git a/src/arome/turb/modd_dynn.F90 b/src/arome/turb/modd_dynn.F90
new file mode 100644
index 0000000000000000000000000000000000000000..f74faf12893795d671d44a546aacf3d3f5a47681
--- /dev/null
+++ b/src/arome/turb/modd_dynn.F90
@@ -0,0 +1,4 @@
+MODULE MODD_DYN_n
+IMPLICIT NONE
+LOGICAL, PARAMETER :: LOCEAN=.FALSE.
+ENDMODULE MODD_DYN_n
diff --git a/src/arome/turb/modd_turbn.F90 b/src/arome/turb/modd_turbn.F90
new file mode 100644
index 0000000000000000000000000000000000000000..3d11a7b591941538033c61d4e9fa3e1f71b1d8d1
--- /dev/null
+++ b/src/arome/turb/modd_turbn.F90
@@ -0,0 +1,3 @@
+MODULE MODD_TURB_n
+ CHARACTER (LEN=4), SAVE  :: CTURBLEN='BL89'
+ENDMODULE MODD_TURB_n
diff --git a/src/arome/turb/modi_bl89.F90 b/src/arome/turb/modi_bl89.F90
index cc334a48540f6cb76f9d2c35a3c9b203dbc87fd6..e451772994675a8640f6cce983666f77da09a8aa 100644
--- a/src/arome/turb/modi_bl89.F90
+++ b/src/arome/turb/modi_bl89.F90
@@ -2,7 +2,7 @@
       MODULE MODI_BL89
 !     ################
 INTERFACE
-      SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PLM)
+      SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
 !
 INTEGER,                  INTENT(IN)  :: KKA 
 INTEGER,                  INTENT(IN)  :: KKU 
@@ -14,6 +14,7 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTHLM
 INTEGER,                  INTENT(IN)  :: KRR
 REAL, DIMENSION(:,:,:,:), INTENT(IN)  :: PRM
 REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTKEM
+REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PSHEAR
 REAL, DIMENSION(:,:,:),   INTENT(OUT) :: PLM
 
 END SUBROUTINE BL89
diff --git a/src/arome/turb/modi_compute_bl89_ml.F90 b/src/arome/turb/modi_compute_bl89_ml.F90
index 7d9c50e82b2076526ba33a1e2677c4e83d5d57d8..c42759094d59acf316dad60bd0f7ae0b13fa7a31 100644
--- a/src/arome/turb/modi_compute_bl89_ml.F90
+++ b/src/arome/turb/modi_compute_bl89_ml.F90
@@ -6,7 +6,7 @@ INTERFACE
 
 !     ###################################################################
       SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
-             PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
+             PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
 !     ###################################################################
 
 !*               1.1  Declaration of Arguments
@@ -16,14 +16,16 @@ INTEGER,                INTENT(IN)   :: KKB          ! near ground physical inde
 INTEGER,                INTENT(IN)   :: KKE          ! uppest atmosphere physical index
 INTEGER,                INTENT(IN)   :: KKU          ! uppest atmosphere array index
 INTEGER,                INTENT(IN)   :: KKL          ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(:,:),   INTENT(IN)  :: PDZZ2D
-REAL, DIMENSION(:),     INTENT(IN)  :: PTKEM_DEP
-REAL, DIMENSION(:),     INTENT(IN)  :: PG_O_THVREF
-REAL, DIMENSION(:,:),   INTENT(IN)  :: PVPT
-INTEGER,                INTENT(IN)  :: KK
-LOGICAL,                INTENT(IN)  :: OUPORDN
+REAL, DIMENSION(:,:),   INTENT(IN)  :: PDZZ2D        ! height difference between two mass levels
+REAL, DIMENSION(:),     INTENT(IN)  :: PTKEM_DEP     ! TKE to consume
+REAL, DIMENSION(:),     INTENT(IN)  :: PG_O_THVREF   ! g/ThetaVRef at the departure point
+REAL, DIMENSION(:,:),   INTENT(IN)  :: PVPT          ! ThetaV on mass levels
+INTEGER,                INTENT(IN)  :: KK            ! index of departure level
+LOGICAL,                INTENT(IN)  :: OUPORDN       ! switch to compute upward (true) or
+                                                     !   downward (false) mixing length
 LOGICAL,                INTENT(IN)  :: OFLUX         ! Computation must be done from flux level
-REAL, DIMENSION(:),     INTENT(OUT) :: PLWORK
+REAL, DIMENSION(:),     INTENT(OUT) :: PLWORK        ! Resulting mixing length
+REAL, DIMENSION(:,:),   INTENT(IN)  :: PSHEAR        ! vertical wind shear for RM17 mixing length
 
 END SUBROUTINE COMPUTE_BL89_ML
 
diff --git a/src/arome/turb/turb.F90 b/src/arome/turb/turb.F90
index 3e52f07833d187f8b94318ebd68269f9b2c15d2a..6b219e33fd20d346984e1538fedd9f8f3df925d5 100644
--- a/src/arome/turb/turb.F90
+++ b/src/arome/turb/turb.F90
@@ -446,6 +446,7 @@ REAL                :: ZALPHA       ! proportionnality constant between Dz/2 and
 !                                   ! BL89 mixing length near the surface
 !
 REAL :: ZTIME1, ZTIME2
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3))::  ZSHEAR, ZDUDZ, ZDVDZ
 !
 !*      1.PRELIMINARIES
 !         -------------
@@ -607,7 +608,8 @@ SELECT CASE (HTURBLEN)
 !           ------------------
 
   CASE ('BL89')
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZLM)
+    ZSHEAR=0.
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZSHEAR,ZLM)
 !
 !*      3.2 Delta mixing length
 !           -------------------
@@ -1516,7 +1518,8 @@ ELSE
 !*         3.1 BL89 mixing length
 !           ------------------
   CASE ('BL89')
-    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZLM_CLOUD)
+    ZSHEAR=0.
+    CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKEM,ZSHEAR,ZLM_CLOUD)
 !
 !*         3.2 Delta mixing length
 !           -------------------
diff --git a/src/common/turb/bl89.f90 b/src/common/turb/bl89.F90
similarity index 96%
rename from src/common/turb/bl89.f90
rename to src/common/turb/bl89.F90
index 303d92d262407a571fbe894e1da0274d1626abf3..57056e4815f6fb65ae1ff202472d58ed1f7810c6 100644
--- a/src/common/turb/bl89.f90
+++ b/src/common/turb/bl89.F90
@@ -136,6 +136,7 @@ IF (CPROGRAM=='AROME ') THEN
     ZZZ    (:,JK)   = PZZ    (:,1,JK)
     ZDZZ   (:,JK)   = PDZZ   (:,1,JK)
     ZTHM   (:,JK)   = PTHLM  (:,1,JK)
+    ZSHEAR (:,JK)   = PSHEAR (:,1,JK)
     ZTKEM  (:,JK)   = PTKEM  (:,1,JK)
     ZG_O_THVREF(:,JK)   = XG/PTHVREF(:,1,JK)
   END DO
@@ -236,7 +237,11 @@ DO JK=IKTB,IKTE
           + sqrt(abs( (XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
           + ( -ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - ZVPT(J1D,JK)) ))**2.0 + &
           2. * ZINTE(J1D) * &
+#ifdef REPRO48
+          ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK)/ ZDZZ(J1D,JKK)))) / &
+#else
           (ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK)/ ZDZZ(J1D,JKK))))) / &
+#endif
           (ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK))
         ZLWORK(J1D)=ZLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2)
         ZINTE(J1D) = ZINTE(J1D) - ZPOTE
@@ -263,6 +268,7 @@ DO JK=IKTB,IKTE
     IF(ZTESTM > 0.) THEN
       ZTESTM=0.
       DO J1D=1,IIU*IJU
+        ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
         !--------- SHEAR + STABILITY -----------
         ZPOTE = ZTEST0* &
                 (ZG_O_THVREF(J1D,JK)*(ZHLVPT(J1D,JKK)-ZVPT(J1D,JK)) &
@@ -278,7 +284,11 @@ DO JK=IKTB,IKTE
           (XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK)   &
             + ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK))) )**2    &
             + 2. * ZINTE(J1D) * &
-            ( ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK))))) / &
+#ifdef REPRO48
+             ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK)))) / &
+#else
+            (ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK))))) / &
+#endif
             (ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK))
         ZLWORK(J1D)=ZLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2)
         ZINTE(J1D) = ZINTE(J1D) - ZPOTE
@@ -294,8 +304,13 @@ DO JK=IKTB,IKTE
     ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10_MNHREAL)
     ZLWORK2=MAX(ZLWORK(J1D),1.E-10_MNHREAL)
     ZPOTE = ZLWORK1 / ZLWORK2
+#ifdef REPRO48
+    ZLWORK2=1.d0 + ZPOTE**(2./3.)
+    ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
+#else
     ZLWORK2=1.d0 + ZPOTE**ZBL89EXP
     ZLM(J1D,JK) = ZLWORK1*(2./ZLWORK2)**ZUSRBL89
+#endif
   END DO
 
 ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
diff --git a/src/common/turb/compute_bl89_ml.f90 b/src/common/turb/compute_bl89_ml.F90
similarity index 97%
rename from src/common/turb/compute_bl89_ml.f90
rename to src/common/turb/compute_bl89_ml.F90
index 5a75011c8ed80ee7d7a3b364136a9ecff36168ea..897bfc12fab3710b1e672d573cc197fa7d24007e 100644
--- a/src/common/turb/compute_bl89_ml.f90
+++ b/src/common/turb/compute_bl89_ml.F90
@@ -127,9 +127,9 @@ IF (OUPORDN.EQV..TRUE.) THEN
      ! Lenght travelled by parcel to nullify energy
      ZLWORK2(J1D)=        ( - PG_O_THVREF(J1D) *                     &
             (  ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) )                              &
-            - XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))) &
+            - XRM17*PSHEAR(J1D,KK)*SQRT(ABS(PTKEM_DEP(J1D))) &
           + SQRT (ABS(                                                       &
-            (XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))) +  &
+            (XRM17*PSHEAR(J1D,KK)*SQRT(ABS(PTKEM_DEP(J1D))) +  &
              PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2  &
             + 2. * ZINTE(J1D) * PG_O_THVREF(J1D)                        &
                  * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) ))    ) /             &
@@ -194,7 +194,7 @@ IF (OUPORDN.EQV..FALSE.) THEN
         ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
         ZTESTM=ZTESTM+ZTEST0
         ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
-      ZLWORK2(J1D)=        ( + PG_O_THVREF(J1D) *                     &
+        ZLWORK2(J1D)=        ( + PG_O_THVREF(J1D) *                     &
             (  PVPT(J1D,JKK) - PVPT(J1D,KK) )                              &
              -XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
           + SQRT (ABS(                                                       &
diff --git a/src/common/turb/rmc01.f90 b/src/common/turb/rmc01.F90
similarity index 100%
rename from src/common/turb/rmc01.f90
rename to src/common/turb/rmc01.F90
diff --git a/src/common/turb/sbl_depth.f90 b/src/common/turb/sbl_depth.F90
similarity index 100%
rename from src/common/turb/sbl_depth.f90
rename to src/common/turb/sbl_depth.F90