diff --git a/src/arome/turb/emoist.F90 b/src/common/turb/emoist.F90
similarity index 87%
rename from src/arome/turb/emoist.F90
rename to src/common/turb/emoist.F90
index ff243f221f02f8d855d0272f67ba426d4c89f079..f17bc3168c165b43e3bac222dad961b974e4834f 100644
--- a/src/arome/turb/emoist.F90
+++ b/src/common/turb/emoist.F90
@@ -1,3 +1,7 @@
+!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.
 !     ######spl
 FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
 USE PARKIND1, ONLY : JPRB
@@ -45,12 +49,14 @@ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !!       J. Stein       Feb  28, 1996   optimization + Doctorization
 !!       J. Stein       Spet 15, 1996   Amoist previously computed
 !!       J.-P. Pinty    May  20, 2003   Improve EMOIST expression
+!!                  03/2021 (JL Redelsperger) Ocean model case 
 !! 
 !! ----------------------------------------------------------------------
 !
 !*       0. DECLARATIONS
 !           ------------
 USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
 !
 IMPLICIT NONE
 !
@@ -84,15 +90,24 @@ INTEGER                               :: JRR     ! moist loop counter
 !*       1. COMPUTE EMOIST
 !           --------------
 !
-!
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('EMOIST',0,ZHOOK_HANDLE)
-IF ( KRR == 0 ) THEN                                ! dry case
-  PEMOIST(:,:,:) = 0.
-ELSE IF ( KRR == 1 ) THEN                           ! only vapor
+!
+IF (LOCEAN) THEN
+ IF ( KRR == 0 ) THEN                                ! Unsalted
+   PEMOIST(:,:,:) = 0.
+ ELSE
+   PEMOIST(:,:,:) = 1.                              ! Salted case
+ END IF
+!
+ELSE
+!
+ IF ( KRR == 0 ) THEN                                ! dry case
+   PEMOIST(:,:,:) = 0.
+ ELSE IF ( KRR == 1 ) THEN                           ! only vapor
   ZDELTA = (XRV/XRD) - 1.
   PEMOIST(:,:,:) = ZDELTA*PTHLM(:,:,:)
-ELSE                                                ! liquid water & ice present
+ ELSE                                                ! liquid water & ice present
   ZDELTA = (XRV/XRD) - 1.
   ZRW(:,:,:) = PRM(:,:,:,1)
 !
@@ -140,8 +155,9 @@ ELSE                                                ! liquid water & ice present
                             / (1. + ZRW(:,:,:))                                &
          ) * PAMOIST(:,:,:) * 2. * PSRCM(:,:,:)
   END IF
-END IF
+ END IF
 !
+END IF
 !---------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('EMOIST',1,ZHOOK_HANDLE)
diff --git a/src/arome/turb/etheta.F90 b/src/common/turb/etheta.F90
similarity index 88%
rename from src/arome/turb/etheta.F90
rename to src/common/turb/etheta.F90
index 77f27da4ef8a5da82a6e977c20a686ba2b82db51..f0506bd89dc5fc455e55281835ba7c1ab6b2365c 100644
--- a/src/arome/turb/etheta.F90
+++ b/src/common/turb/etheta.F90
@@ -1,3 +1,7 @@
+!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.
 !     ######spl
 FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
 USE PARKIND1, ONLY : JPRB
@@ -45,12 +49,13 @@ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
 !!       J. Stein       Feb  28, 1996   optimization + Doctorization
 !!       J. Stein       Sept 15, 1996   Atheta previously computed
 !!       J.-P. Pinty    May  20, 2003   Improve ETHETA expression
-!!
+!!       J.L Redelsperger    03, 2021   Ocean Model Case 
 !! ----------------------------------------------------------------------
 !
 !*       0. DECLARATIONS
 !           ------------
 USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
 !
 IMPLICIT NONE
 !
@@ -90,12 +95,15 @@ INTEGER                               :: JRR     ! moist loop counter
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('ETHETA',0,ZHOOK_HANDLE)
-IF ( KRR == 0 ) THEN                                ! dry case
+IF (LOCEAN) THEN                                    ! ocean case
+   PETHETA(:,:,:) =  1.
+ELSE   
+ IF ( KRR == 0) THEN                                ! dry case
   PETHETA(:,:,:) = 1.
-ELSE IF ( KRR == 1 ) THEN                           ! only vapor
+ ELSE IF ( KRR == 1 ) THEN                           ! only vapor
   ZDELTA = (XRV/XRD) - 1.
   PETHETA(:,:,:) = 1. + ZDELTA*PRM(:,:,:,1)
-ELSE                                                ! liquid water & ice present
+ ELSE                                                ! liquid water & ice present
   ZDELTA = (XRV/XRD) - 1.
   ZRW(:,:,:) = PRM(:,:,:,1)
 !
@@ -138,8 +146,9 @@ ELSE                                                ! liquid water & ice present
                             / (1. + ZRW(:,:,:))                                &
          ) * PATHETA(:,:,:) * 2. * PSRCM(:,:,:)
   END IF
-END IF
+ END IF
 !
+END IF
 !---------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('ETHETA',1,ZHOOK_HANDLE)
diff --git a/src/arome/turb/shuman_mf.F90 b/src/common/turb/shuman_mf.F90
similarity index 98%
rename from src/arome/turb/shuman_mf.F90
rename to src/common/turb/shuman_mf.F90
index 56d09f4428201d4d11e4f8790450bfaa05576c7f..4720135892d6c4c6b0a1122a76d99287ba2bbc44 100644
--- a/src/arome/turb/shuman_mf.F90
+++ b/src/common/turb/shuman_mf.F90
@@ -1,3 +1,7 @@
+!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.
 !     ##################
       MODULE MODI_SHUMAN_MF
 !     ##################
diff --git a/src/mesonh/turb/emoist.f90 b/src/mesonh/turb/emoist.f90
deleted file mode 100644
index 7703fb388efc28be14bffeb7d070dcdc0167af19..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/emoist.f90
+++ /dev/null
@@ -1,185 +0,0 @@
-!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 MODI_EMOIST
-!#################
-!
-INTERFACE
-!
-FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
-!
-INTEGER                              :: KRR        ! number of moist var.
-INTEGER                              :: KRRI       ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM    ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM      ! Mixing ratios, where
-!                                    PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PAMOIST   ! Amoist
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PEMOIST ! result
-!
-END FUNCTION EMOIST
-!
-END INTERFACE
-!
-END MODULE MODI_EMOIST
-!
-!   ############################################################################
-FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
-!   ############################################################################
-!
-!      PURPOSE
-!!     -------
-!      EMOIST computes the coefficient Emoist in the flottability turbulent
-!      flux. This coefficient relates the vertical flux of the virtual potential
-!      temperature ( <Thv' W'> ) to the vertical flux of the conservative mixing
-!      ratio ( <Rnp' W'> ). 
-!
-!!**   METHOD
-!!     ------
-!!     The virtual potential temperature perturbation is linearized in function
-!!     of Thl' and Rnp'. The result is
-!!        Thv'= ( ZA + ZC * Atheta * 2 * SRC ) Thl' 
-!!             +( ZB + ZC * Amoist * 2 * SRC ) Rnp'
-!!     From this relation, we can compute the verical turbulent fluxes.
-!! 
-!!     EXTERNAL
-!!     --------
-!!       NONE
-!!
-!!     IMPLICIT ARGUMENTS
-!!     ------------------
-!!       Module MODD_CST : contains physical constants.
-!!         XRV, XRD  : R for water vapor and dry air
-!!   
-!!     REFERENCE
-!!     ---------
-!! 
-!!       NONE
-!!
-!!
-!!     AUTHOR
-!!     ------
-!!       Jean-Marie Carriere      * Meteo-France *
-!!
-!!     MODIFICATIONS
-!!     -------------
-!!       Original       20/03/95
-!!
-!!       J. Stein       Feb  28, 1996   optimization + Doctorization
-!!       J. Stein       Spet 15, 1996   Amoist previously computed
-!!       J.-P. Pinty    May  20, 2003   Improve EMOIST expression
-!!                  03/2021 (JL Redelsperger) Ocean model case 
-!! 
-!! ----------------------------------------------------------------------
-!
-!*       0. DECLARATIONS
-!           ------------
-USE MODD_CST
-USE MODD_DYN_n, ONLY : LOCEAN
-!
-IMPLICIT NONE
-!
-!*       0.1 declarations of arguments and result
-!
-!
-INTEGER                              :: KRR        ! number of moist var.
-INTEGER                              :: KRRI       ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM    ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM      ! Mixing ratios, where
-!                                    PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PAMOIST   ! Amoist
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PEMOIST ! result
-!
-!*       0.2 declarations of local variables
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::       &
-                                        ZA, ZRW
-!                ZA = coeft A, ZRW = total mixing ratio rw 
-REAL                                  :: ZDELTA  ! = Rv/Rd - 1
-INTEGER                               :: JRR     ! moist loop counter
-!
-!---------------------------------------------------------------------------
-!
-!
-!*       1. COMPUTE EMOIST
-!           --------------
-IF (LOCEAN) THEN
- IF ( KRR == 0 ) THEN                                ! Unsalted
-   PEMOIST(:,:,:) = 0.
- ELSE
-   PEMOIST(:,:,:) = 1.                              ! Salted case
- END IF
-!
-ELSE
-!
- IF ( KRR == 0 ) THEN                                ! dry case
-   PEMOIST(:,:,:) = 0.
- ELSE IF ( KRR == 1 ) THEN                           ! only vapor
-  ZDELTA = (XRV/XRD) - 1.
-  PEMOIST(:,:,:) = ZDELTA*PTHLM(:,:,:)
- ELSE                                                ! liquid water & ice present
-  ZDELTA = (XRV/XRD) - 1.
-  ZRW(:,:,:) = PRM(:,:,:,1)
-!
-  IF ( KRRI>0) THEN  ! rc and ri case  
-    ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,3)
-    DO JRR=5,KRR
-      ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,JRR)
-    ENDDO
-    ZA(:,:,:) = 1. + (                                    &  ! Compute A
-              (1.+ZDELTA) * (PRM(:,:,:,1) - PRM(:,:,:,2) - PRM(:,:,:,4)) &
-              -ZRW(:,:,:)                                                &
-                     )  /  (1. + ZRW(:,:,:)) 
-  !
-  !   Emoist = ZB + ZC * Amoist
-  !   ZB is computed from line 1 to line 2
-  !   ZC is computed from line 3 to line 5
-  !   Amoist* 2 * SRC is computed at line 6
-  !
-    PEMOIST(:,:,:) = ZDELTA * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*(               &
-                                                    PRM(:,:,:,2)+PRM(:,:,:,4)))&
-                            / (1. + ZRW(:,:,:))                                &
-        +( PLOCPEXNM(:,:,:) * ZA(:,:,:)                                        &
-               -(1.+ZDELTA) * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*(               &
-                                                    PRM(:,:,:,2)+PRM(:,:,:,4)))&
-                            / (1. + ZRW(:,:,:))                                &
-         ) * PAMOIST(:,:,:) * 2. * PSRCM(:,:,:)
-  ELSE
-    DO JRR=3,KRR
-      ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,JRR)
-    ENDDO
-    ZA(:,:,:) = 1. + (                                    &  ! Compute ZA
-              (1.+ZDELTA) * (PRM(:,:,:,1) - PRM(:,:,:,2)) &
-              -ZRW(:,:,:)                                 &
-                     )  /  (1. + ZRW(:,:,:)) 
-  !
-  !   Emoist = ZB + ZC * Amoist
-  !   ZB is computed from line 1 to line 2
-  !   ZC is computed from line 3 to line 5
-  !   Amoist* 2 * SRC is computed at line 6
-  !
-    PEMOIST(:,:,:) = ZDELTA * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*PRM(:,:,:,2))   &
-                            / (1. + ZRW(:,:,:))                                &
-        +( PLOCPEXNM(:,:,:) * ZA(:,:,:)                                        &
-               -(1.+ZDELTA) * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*PRM(:,:,:,2))   &
-                            / (1. + ZRW(:,:,:))                                &
-         ) * PAMOIST(:,:,:) * 2. * PSRCM(:,:,:)
-  END IF
- END IF
-!
-END IF
-!---------------------------------------------------------------------------
-!
-END FUNCTION EMOIST
diff --git a/src/mesonh/turb/etheta.f90 b/src/mesonh/turb/etheta.f90
deleted file mode 100644
index 3ef29178b721660ec33639f4199166e0f6d1a9d0..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/etheta.f90
+++ /dev/null
@@ -1,180 +0,0 @@
-!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 MODI_ETHETA
-!#################
-!
-INTERFACE
-!
-FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
-!
-INTEGER                              :: KRR          ! number of moist var.
-INTEGER                              :: KRRI         ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM     ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM       ! Mixing ratios, where
-!                                      PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PATHETA   ! Atheta
-!                                                    
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PETHETA ! result
-!
-!
-END FUNCTION ETHETA
-!
-END INTERFACE
-!
-END MODULE MODI_ETHETA
-!
-!   ############################################################################
-FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
-!   ############################################################################
-!
-!      PURPOSE
-!!     -------
-!      ETHETA computes the coefficient Etheta in the flottability turbulent
-!      flux. This coefficient relates the vertical flux of the virtual potential
-!      temperature ( <Thv' W'> ) to the vertical flux of the conservative potential
-!      temperature ( <Thl' W'> ).
-!
-!!**   METHOD
-!!     ------
-!!
-!!     The virtual potential temperature perturbation is linearized in function
-!!     of Thl' and Rnp'. The result is
-!!        Thv'= ( ZA + ZC * Atheta * 2 * SRC ) Thl' 
-!!             +( ZB + ZC * Amoist * 2 * SRC ) Rnp'
-!!     From this relation, we can compute the vertical turbulent fluxes.
-!!
-!!     EXTERNAL
-!!     --------
-!!
-!!        NONE
-!!
-!!     IMPLICIT ARGUMENTS
-!!     ------------------
-!!       Module MODD_CST : contains physical constants.
-!!         XRV, XRD  : R for water vapor and dry air
-!!            
-!!     REFERENCE
-!!     ---------
-!!
-!!
-!!     AUTHOR
-!!     ------
-!!       Jean-Marie Carriere      * Meteo-France *
-!!
-!!     MODIFICATIONS
-!!     -------------
-!!       Original       20/03/95
-!!     
-!!       J. Stein       Feb  28, 1996   optimization + Doctorization
-!!       J. Stein       Sept 15, 1996   Atheta previously computed
-!!       J.-P. Pinty    May  20, 2003   Improve ETHETA expression
-!!       J.L Redelsperger    03, 2021   Ocean Model Case 
-!! ----------------------------------------------------------------------
-!
-!*       0. DECLARATIONS
-!           ------------
-USE MODD_CST
-USE MODD_DYN_n, ONLY : LOCEAN
-!
-IMPLICIT NONE
-!
-!*       0.1 declarations of arguments and result
-!
-!
-INTEGER                              :: KRR          ! number of moist var.
-INTEGER                              :: KRRI         ! number of ice var.
-!
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PTHLM     ! Conservative pot. temperature
-REAL, DIMENSION(:,:,:,:), INTENT(IN) ::   PRM       ! Mixing ratios, where
-!                                      PRM(:,:,:,1) = conservative mixing ratio
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PLOCPEXNM ! Lv(T)/Cp/Exner at time t-1
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PATHETA   ! Atheta
-!                                                    
-REAL, DIMENSION(:,:,:),  INTENT(IN)  ::   PSRCM     ! Normalized 2dn_order
-                                                    ! moment s'r'c/2Sigma_s2
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)):: PETHETA ! result
-!
-!
-!
-!*       0.2 declarations of local variables
-!
-REAL,DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::       &
-                                        ZA, ZRW
-!                ZA = coeft A, ZRW = total mixing ratio rw
-REAL                                  :: ZDELTA  ! = Rv/Rd - 1
-INTEGER                               :: JRR     ! moist loop counter
-!
-!---------------------------------------------------------------------------
-!
-!
-!*       1. COMPUTE ETHETA
-!           --------------
-!
-!
-IF (LOCEAN) THEN                                    ! ocean case
-   PETHETA(:,:,:) =  1.
-ELSE   
- IF ( KRR == 0.) THEN                                ! dry case
-  PETHETA(:,:,:) = 1.
- ELSE IF ( KRR == 1 ) THEN                           ! only vapor
-  ZDELTA = (XRV/XRD) - 1.
-  PETHETA(:,:,:) = 1. + ZDELTA*PRM(:,:,:,1)
- ELSE                                                ! liquid water & ice present
-  ZDELTA = (XRV/XRD) - 1.
-  ZRW(:,:,:) = PRM(:,:,:,1)
-!
-  IF ( KRRI>0 ) THEN  ! rc and ri case
-    ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,3)
-    DO JRR=5,KRR
-      ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,JRR)
-    ENDDO
-    ZA(:,:,:) = 1. + (                                    &  ! Compute A
-              (1.+ZDELTA) * (PRM(:,:,:,1) - PRM(:,:,:,2) - PRM(:,:,:,4)) &
-              -ZRW(:,:,:)                                                &
-                     )  /  (1. + ZRW(:,:,:))
-  !
-  !   Etheta = ZA + ZC * Atheta  
-  !   ZC is computed from line 2 to line 5
-  !   - Atheta * 2. * SRC is computed at line 6 
-  !
-    PETHETA(:,:,:) = ZA(:,:,:)                                                 &
-        +( PLOCPEXNM(:,:,:) * ZA(:,:,:)                                        &
-               -(1.+ZDELTA) * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*(               &
-                                                    PRM(:,:,:,2)+PRM(:,:,:,4)))&
-                            / (1. + ZRW(:,:,:))                                &
-         ) * PATHETA(:,:,:) * 2. * PSRCM(:,:,:)
-  ELSE
-    DO JRR=3,KRR
-      ZRW(:,:,:) = ZRW(:,:,:) + PRM(:,:,:,JRR)
-    ENDDO
-    ZA(:,:,:) = 1. + (                                    &  ! Compute A
-              (1.+ZDELTA) * (PRM(:,:,:,1) - PRM(:,:,:,2)) &
-              -ZRW(:,:,:)                                 &
-                     )  /  (1. + ZRW(:,:,:))
-  !
-  !   Etheta = ZA + ZC * Atheta  
-  !   ZC is computed from line 2 to line 5
-  !   - Atheta * 2. * SRC is computed at line 6 
-  !
-    PETHETA(:,:,:) = ZA(:,:,:)                                                 &
-        +( PLOCPEXNM(:,:,:) * ZA(:,:,:)                                        &
-               -(1.+ZDELTA) * (PTHLM(:,:,:) + PLOCPEXNM(:,:,:)*PRM(:,:,:,2))   &
-                            / (1. + ZRW(:,:,:))                                &
-         ) * PATHETA(:,:,:) * 2. * PSRCM(:,:,:)
-  END IF
- END IF
-!
-END IF
-!---------------------------------------------------------------------------
-!
-END FUNCTION ETHETA
diff --git a/src/mesonh/turb/shuman_mf.f90 b/src/mesonh/turb/shuman_mf.f90
deleted file mode 100644
index ce9cde0512e4e7c8ebec82d323448dc1dd2e3c6e..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/shuman_mf.f90
+++ /dev/null
@@ -1,445 +0,0 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
-!MNH_LIC for details. version 1.
-!     ##################
-      MODULE MODI_SHUMAN_MF
-!     ##################
-!
-INTERFACE
-!
-FUNCTION DZF_MF(KKA,KKU,KKL,PA)  RESULT(PDZF)
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux
-                                                 !  side
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZF   ! result at mass
-                                                 ! localization 
-END FUNCTION DZF_MF
-!
-FUNCTION DZM_MF(KKA,KKU,KKL,PA)  RESULT(PDZM)
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass
-                                                 ! localization
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZM   ! result at flux
-                                                 ! side
-END FUNCTION DZM_MF
-!
-FUNCTION MZF_MF(KKA,KKU,KKL,PA)  RESULT(PMZF)
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux
-                                                 !  side
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZF   ! result at mass
-                                                 ! localization 
-END FUNCTION MZF_MF
-!
-FUNCTION MZM_MF(KKA,KKU,KKL,PA)  RESULT(PMZM)
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZM   ! result at flux localization 
-END FUNCTION MZM_MF
-!
-FUNCTION GZ_M_W_MF(KKA,KKU,KKL,PY,PDZZ) RESULT(PGZ_M_W)
-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
-REAL, DIMENSION(:,:), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
-REAL, DIMENSION(:,:), INTENT(IN)  :: PY   ! variable at mass localization
-REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2)) :: PGZ_M_W  ! result at flux side
-END FUNCTION GZ_M_W_MF
-!
-END INTERFACE
-!
-END MODULE MODI_SHUMAN_MF
-!
-!     ###############################
-      FUNCTION MZF_MF(KKA,KKU,KKL,PA)  RESULT(PMZF)
-!     ###############################
-!
-!!****  *MZF* -  SHUMAN_MF operator : mean operator in z direction for a 
-!!                                 variable at a flux side
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this function  is to compute a mean 
-!     along the z direction (K index) for a field PA localized at a z-flux
-!     point (w point). The result is localized at a mass point.
-!
-!!**  METHOD
-!!    ------ 
-!!        The result PMZF(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k+1))
-!!        At k=size(PA,3), PMZF(:,:,k) is defined by PA(:,:,k).
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      NONE
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation of Meso-NH (SHUMAN_MF operators)
-!!      Technical specifications Report of The Meso-NH (chapters 3)
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!	V. Ducrocq       * Meteo France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    04/07/94 
-!!                   optimisation                 20/08/00 J. Escobar
-!!      S. Riette, Jan 2012: Simplification and suppression of array overflow
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of argument and result
-!              ------------------------------------
-!
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux
-                                                 !  side
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZF   ! result at mass
-                                                 ! localization 
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-INTEGER :: JK             ! Loop index in z direction
-!            
-!
-!-------------------------------------------------------------------------------
-!
-!*       1.    DEFINITION OF MZF
-!              ------------------
-!
-DO JK=2,SIZE(PA,2)-1
-  PMZF(:,JK) = 0.5*( PA(:,JK)+PA(:,JK+KKL) )
-END DO
-PMZF(:,KKA) = 0.5*( PA(:,KKA)+PA(:,KKA+KKL) )
-PMZF(:,KKU) = PA(:,KKU)
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION MZF_MF
-!     ###############################
-      FUNCTION MZM_MF(KKA,KKU,KKL,PA)  RESULT(PMZM)
-!     ###############################
-!
-!!****  *MZM* -  SHUMAN_MF operator : mean operator in z direction for a 
-!!                                 mass variable 
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this function  is to compute a mean
-!     along the z direction (K index) for a field PA localized at a mass
-!     point. The result is localized at a z-flux point (w point).
-!
-!!**  METHOD
-!!    ------ 
-!!        The result PMZM(:,:,k) is defined by 0.5*(PA(:,:,k)+PA(:,:,k-1))
-!!        At k=1, PMZM(:,:,1) is defined by PA(:,:,1).
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      NONE
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation of Meso-NH (SHUMAN_MF operators)
-!!      Technical specifications Report of The Meso-NH (chapters 3)  
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!	V. Ducrocq       * Meteo France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    04/07/94 
-!!                   optimisation                 20/08/00 J. Escobar
-!!      S. Riette, Jan 2012: Simplification and suppression of array overflow
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of argument and result
-!              ------------------------------------
-!
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass localization
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PMZM   ! result at flux localization 
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-INTEGER :: JK             ! Loop index in z direction
-!           
-!
-!-------------------------------------------------------------------------------
-!
-!*       1.    DEFINITION OF MZM
-!              ------------------
-!
-DO JK=2,SIZE(PA,2)-1
-  PMZM(:,JK) = 0.5*( PA(:,JK)+PA(:,JK-KKL) )
-END DO
-PMZM(:,KKA) = PA(:,KKA)
-PMZM(:,KKU) = 0.5*( PA(:,KKU)+PA(:,KKU-KKL) )
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION MZM_MF
-!     ###############################
-      FUNCTION DZF_MF(KKA,KKU,KKL,PA)  RESULT(PDZF)
-!     ###############################
-!
-!!****  *DZF* -  SHUMAN_MF operator : finite difference operator in z direction
-!!                                  for a variable at a flux side
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this function  is to compute a finite difference 
-!     along the z direction (K index) for a field PA localized at a z-flux
-!     point (w point). The result is localized at a mass point.
-!
-!!**  METHOD
-!!    ------ 
-!!        The result PDZF(:,:,k) is defined by (PA(:,:,k+1)-PA(:,:,k))
-!!        At k=size(PA,3), PDZF(:,:,k) is defined by 0.
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      NONE
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation of Meso-NH (SHUMAN_MF operators)
-!!      Technical specifications Report of The Meso-NH (chapters 3)  
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!	V. Ducrocq       * Meteo France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    05/07/94 
-!!                   optimisation                 20/08/00 J. Escobar
-!!      S. Riette, Jan 2012: Simplification and suppression of array overflow
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of argument and result
-!              ------------------------------------
-!
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at flux
-                                                 !  side
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZF   ! result at mass
-                                                 ! localization 
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-INTEGER :: JK           ! Loop index in z direction
-!         
-!-------------------------------------------------------------------------------
-!
-!*       1.    DEFINITION OF DZF
-!              ------------------
-!
-DO JK=2,SIZE(PA,2)-1
-  PDZF(:,JK) = PA(:,JK+KKL) - PA(:,JK)
-END DO
-PDZF(:,KKA) = PA(:,KKA+KKL) - PA(:,KKA)
-PDZF(:,KKU) = 0.
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION DZF_MF
-!     ###############################
-      FUNCTION DZM_MF(KKA,KKU,KKL,PA)  RESULT(PDZM)
-!     ###############################
-!
-!!****  *DZM* -  SHUMAN_MF operator : finite difference operator in z direction
-!!                                  for a variable at a mass localization
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this function  is to compute a finite difference 
-!     along the z direction (K index) for a field PA localized at a mass
-!     point. The result is localized at a z-flux point (w point).
-!
-!!**  METHOD
-!!    ------ 
-!!        The result PDZM(:,j,:) is defined by (PA(:,:,k)-PA(:,:,k-1))
-!!        At k=1, PDZM(:,:,k) is defined by 0.
-!!    
-!!
-!!    EXTERNAL
-!!    --------
-!!      NONE
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation of Meso-NH (SHUMAN_MF operators)
-!!      Technical specifications Report of The Meso-NH (chapters 3)  
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!	V. Ducrocq       * Meteo France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    05/07/94 
-!!                   optimisation                 20/08/00 J. Escobar
-!!      S. Riette, Jan 2012: Simplification and suppression of array overflow
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of argument and result
-!              ------------------------------------
-!
-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
-REAL, DIMENSION(:,:), INTENT(IN)       :: PA     ! variable at mass
-                                                 ! localization
-REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: PDZM   ! result at flux
-                                                 ! side
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-INTEGER :: JK            ! Loop index in z direction
-!           
-!-------------------------------------------------------------------------------
-!
-!*       1.    DEFINITION OF DZM
-!              ------------------
-!
-DO JK=2,SIZE(PA,2)-1
-  PDZM(:,JK) = PA(:,JK) - PA(:,JK-KKL)
-END DO
-PDZM(:,KKA) = 0.
-PDZM(:,KKU) = PA(:,KKU) - PA(:,KKU-KKL)
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION DZM_MF
-
-!     ###############################
-      FUNCTION GZ_M_W_MF(KKA,KKU,KKL,PY,PDZZ) RESULT(PGZ_M_W)
-!     ###############################
-!
-!!****  *GZ_M_W * - Compute the gradient along z direction for a
-!!       variable localized at a mass point
-!!
-!!    PURPOSE
-!!    -------
-!
-!!**  METHOD
-!!    ------
-!                    dzm(PY)
-!       PGZ_M_W =    -------
-!                     d*zz
-!!
-!!
-!!    EXTERNAL
-!!    --------
-!!      NONE
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!    S.Riette moving of code previously in compute_mf_cloud code
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    25 Aug 2011
-!!      S. Riette, Jan 2012: Simplification and suppression of array overflow
-!!
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-!!
-!
-!-------------------------------------------------------------------------------
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of argument and result
-!              ------------------------------------
-!
-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
-REAL, DIMENSION(:,:), INTENT(IN)  :: PDZZ ! Metric coefficient d*zz
-REAL, DIMENSION(:,:), INTENT(IN)  :: PY   ! variable at mass localization
-REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2)) :: PGZ_M_W  ! result at flux side
-!
-!*       0.2   Declarations of local variables
-!              -------------------------------
-!
-INTEGER  JK
-!-------------------------------------------------------------------------------
-!
-!*       1.    COMPUTE THE GRADIENT ALONG Z
-!              -----------------------------
-!
-DO JK=2,SIZE(PY,2)-1
-  PGZ_M_W(:,JK) = (PY(:,JK) - PY(:,JK-KKL)) / PDZZ(:,JK)
-END DO
-PGZ_M_W(:,KKA) = 0.
-PGZ_M_W(:,KKU) = (PY(:,KKU) - PY(:,KKU-KKL)) / PDZZ(:,KKU)
-!
-!-------------------------------------------------------------------------------
-!
-END FUNCTION GZ_M_W_MF