diff --git a/src/common/aux/gradient_m_phy.F90 b/src/common/aux/gradient_m_phy.F90
index fb2119e05789f0f43dddb8d05c45213cefde1eb8..94dd54aea0e555e5e0ad8ecaaa5e9c0ef0dc6c13 100644
--- a/src/common/aux/gradient_m_phy.F90
+++ b/src/common/aux/gradient_m_phy.F90
@@ -99,4 +99,251 @@ PGZ_M_W(IIB:IIE,IJB:IJE,D%NKA)= PGZ_M_W(IIB:IIE,IJB:IJE,D%NKU) ! -999.
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE GZ_M_W_PHY
+!
+SUBROUTINE GX_M_M_PHY(D,OFLAT,PA,PDXX,PDZZ,PDZX,PGX_M_M)
+      USE PARKIND1, ONLY : JPRB
+      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!     #######################################################
+!
+!!****  *GX_M_M* - Cartesian Gradient operator: 
+!!                          computes the gradient in the cartesian X
+!!                          direction for a variable placed at the 
+!!                          mass point and the result is placed at
+!!                          the mass point.
+!!    PURPOSE
+!!    -------
+!       The purpose of this function is to compute the discrete gradient 
+!     along the X cartesian direction for a field PA placed at the 
+!     mass point. The result is placed at the mass point.
+!
+!
+!                       (          ______________z )
+!                       (          (___x         ) )
+!                    1  (    _x    (d*zx dzm(PA) ) ) 
+!      PGX_M_M =   ---- (dxf(PA) - (------------)) )
+!                  ___x (          (             ) )
+!                  d*xx (          (      d*zz   ) )     
+!
+!       
+!
+!!**  METHOD
+!!    ------
+!!      The Chain rule of differencing is applied to variables expressed
+!!    in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!!    the cartesian system
+!!        
+!!    EXTERNAL
+!!    --------
+!!      MXM,MXF,MZF     : Shuman functions (mean operators)
+!!      DXF,DZF         : Shuman functions (finite difference operators)
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      MODD_CONF : LFLAT
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!!      A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!!    AUTHOR
+!!    ------
+!!      Joan Cuxart        *INM and Meteo-France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    18/07/94
+!!                  19/07/00  add the LFLAT switch (J. Stein)
+!!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!-------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE SHUMAN_PHY, ONLY: DXF_PHY, MZF_PHY, DZM_PHY, MXF_PHY, MXM_PHY
+!
+IMPLICIT NONE
+!
+!
+!*       0.1   declarations of arguments and result
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PA      ! variable at the mass point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDXX    ! metric coefficient dxx
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDZX    ! metric coefficient dzx
+LOGICAL, INTENT(IN) :: OFLAT
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGX_M_M ! result mass point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4, ZWORK5, ZWORK6, ZMXF_PDXX
+!
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!*       0.2   declaration of local variables
+!
+!              NONE
+!
+!----------------------------------------------------------------------------
+!
+!*       1.    DEFINITION of GX_M_M
+!              --------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GX_M_M',0,ZHOOK_HANDLE)
+!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+CALL MXF_PHY(D,PDXX,ZMXF_PDXX)
+CALL MXM_PHY(D,PA,ZWORK1)
+CALL DXF_PHY(D,ZWORK1,ZWORK2)
+!
+IF (.NOT. OFLAT) THEN
+  CALL DZM_PHY(D,PA,ZWORK3)
+  CALL MXF_PHY(D,PDZX,ZWORK4)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  CALL MZF_PHY(D,ZWORK5,ZWORK6)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  PGX_M_M(IIB:IIE,IJB:IJE,1:D%NKT)= (ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) - ZWORK6(IIB:IIE,IJB:IJE,1:D%NKT)) &
+                                    / ZMXF_PDXX(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+ELSE
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  PGX_M_M(IIB:IIE,IJB:IJE,1:D%NKT)= ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) / ZMXF_PDXX(IIB:IIE,IJB:IJE,1:D%NKT) 
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+END IF
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GX_M_M',1,ZHOOK_HANDLE)
+END SUBROUTINE GX_M_M_PHY
+!
+      SUBROUTINE GY_M_M_PHY(D,OFLAT,PA,PDYY,PDZZ,PDZY,PGY_M_M)
+      USE PARKIND1, ONLY : JPRB
+      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!     #######################################################
+!
+!!****  *GY_M_M* - Cartesian Gradient operator: 
+!!                          computes the gradient in the cartesian Y
+!!                          direction for a variable placed at the 
+!!                          mass point and the result is placed at
+!!                          the mass point.
+!!    PURPOSE
+!!    -------
+!       The purpose of this function is to compute the discrete gradient 
+!     along the Y cartesian direction for a field PA placed at the 
+!     mass point. The result is placed at the mass point.
+!
+!
+!                       (          ______________z )
+!                       (          (___y         ) )
+!                    1  (    _y    (d*zy dzm(PA) ) ) 
+!      PGY_M_M =   ---- (dyf(PA) - (------------)) )
+!                  ___y (          (             ) )
+!                  d*yy (          (      d*zz   ) )     
+!
+!       
+!!**  METHOD
+!!    ------
+!!      The Chain rule of differencing is applied to variables expressed
+!!    in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!!    the cartesian system
+!!        
+!!    EXTERNAL
+!!    --------
+!!      MYM,MYF,MZF     : Shuman functions (mean operators)
+!!      DYF,DZF         : Shuman functions (finite difference operators)
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      MODD_CONF : LFLAT
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!!      A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!!    AUTHOR
+!!    ------
+!!      Joan Cuxart        *INM and Meteo-France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    18/07/94
+!!                  19/07/00  add the LFLAT switch (J. Stein)
+!-------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE SHUMAN_PHY, ONLY: DYF_PHY, MZF_PHY, DZM_PHY, MYF_PHY, MYM_PHY
+!
+IMPLICIT NONE
+!
+!
+!*       0.1   declarations of arguments and result
+!
+TYPE(DIMPHYEX_t),        INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PA      ! variable at the mass point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDYY    ! metric coefficient dyy
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),  INTENT(IN)  :: PDZY    ! metric coefficient dzy
+LOGICAL, INTENT(IN) :: OFLAT
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT),INTENT(OUT) :: PGY_M_M ! result mass point
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4, ZWORK5, ZMYF_PDYY
+!
+INTEGER :: IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!
+!*       0.2   declaration of local variables
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('GY_M_M',0,ZHOOK_HANDLE)
+!
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+!----------------------------------------------------------------------------
+!
+!*       1.    DEFINITION of GY_M_M
+!              --------------------
+!
+CALL MYM_PHY(D,PA,ZWORK1)
+CALL DYF_PHY(D,ZWORK1,ZWORK2)
+CALL MYF_PHY(D,PDYY,ZMYF_PDYY)
+!
+IF (.NOT. OFLAT) THEN
+  !
+  CALL DZM_PHY(D,PA,ZWORK3)
+  CALL MYF_PHY(D,PDZY,ZWORK4)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                   / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  CALL MZF_PHY(D,ZWORK5,ZWORK4)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  PGY_M_M(IIB:IIE,IJB:IJE,1:D%NKT)= (ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)-ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)) &
+                                    /ZMYF_PDYY(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+ELSE
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  PGY_M_M(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ZMYF_PDYY(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+ENDIF  
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('GY_M_M',1,ZHOOK_HANDLE)
+END SUBROUTINE GY_M_M_PHY
 END MODULE MODE_GRADIENT_M_PHY
diff --git a/src/common/aux/shuman_phy.F90 b/src/common/aux/shuman_phy.F90
index a81e6f216d47630fd61afe02c780c4d659b68c85..dbe34cd2c5e75668a485618de77cc2bde574572a 100644
--- a/src/common/aux/shuman_phy.F90
+++ b/src/common/aux/shuman_phy.F90
@@ -156,8 +156,8 @@ IMPLICIT NONE
 !              ------------------------------------
 !
 TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
-REAL, DIMENSION(:,:,:), INTENT(IN)  :: PA     ! variable at mass localization
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMYM   ! result at flux localization 
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMYM   ! result at flux localization 
 !
 !*       0.2   Declarations of local variables
 !              -------------------------------
@@ -411,8 +411,8 @@ IMPLICIT NONE
 !              ------------------------------------
 !
 TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
-REAL, DIMENSION(:,:,:), INTENT(IN)  :: PA     ! variable at mass localization
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMXM   ! result at flux localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)  :: PA     ! variable at mass localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PMXM   ! result at flux localization
 !
 !*       0.2   Declarations of local variables
 !              -------------------------------
@@ -697,4 +697,158 @@ PDZF(:,:,D%NKU)    = -999.
 !
 IF (LHOOK) CALL DR_HOOK('DZF',1,ZHOOK_HANDLE)
 END SUBROUTINE DZF_PHY
+!
+!     ###############################
+      SUBROUTINE DXF_PHY(D,PA,PDXF)
+      USE PARKIND1, ONLY : JPRB
+      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!     ###############################
+!
+!!****  *DXF* -  Shuman operator : finite difference operator in x direction
+!!                                  for a variable at a flux side
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function  is to compute a finite difference
+!     along the x direction (I index) for a field PA localized at a x-flux
+!     point (u point). The result is localized at a mass point.
+!
+!!**  METHOD
+!!    ------
+!!        The result PDXF(i,:,:) is defined by (PA(i+1,:,:)-PA(i,:,:))
+!!        At i=size(PA,1), PDXF(i,:,:) are replaced by the values of PDXF,
+!!    which are the right values in the x-cyclic case
+!!
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_PARAMETERS: declaration of parameter variables
+!!        JPHEXT: define the number of marginal points out of the
+!!        physical domain along the horizontal directions.
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (SHUMAN operators)
+!!      Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!      V. Ducrocq       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/07/94
+!!      Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of argument and result
+!              ------------------------------------
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)  :: PA     ! variable at flux
+                                                          !  side
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PDXF   ! result at mass
+                                                          ! localization
+!
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('DXF',0,ZHOOK_HANDLE)
+!
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('DXF',1,ZHOOK_HANDLE)
+END SUBROUTINE DXF_PHY
+!
+!     ###############################
+      SUBROUTINE DYF_PHY(D,PA,PDYF)
+      USE PARKIND1, ONLY : JPRB
+      USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!     ###############################
+!
+!!****  *DYF* -  Shuman operator : finite difference operator in y direction
+!!                                  for a variable at a flux side
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function  is to compute a finite difference
+!     along the y direction (J index) for a field PA localized at a y-flux
+!     point (v point). The result is localized at a mass point.
+!
+!!**  METHOD
+!!    ------
+!!        The result PDYF(:,j,:) is defined by (PA(:,j+1,:)-PA(:,j,:))
+!!        At j=size(PA,2), PDYF(:,j,:) are replaced by the values of PDYM,
+!!    which are the right values in the y-cyclic case
+!!
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_PARAMETERS: declaration of parameter variables
+!!        JPHEXT: define the number of marginal points out of the
+!!        physical domain along the horizontal directions.
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (SHUMAN operators)
+!!      Technical specifications Report of The Meso-NH (chapters 3)
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!      V. Ducrocq       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/07/94
+!!      Modification to include the periodic case 13/10/94 J.Stein
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of argument and result
+!              ------------------------------------
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)  :: PA     ! variable at flux
+                                                          !  side
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PDYF   ! result at mass
+                                                          ! localization
+!-------------------------------------------------------------------------------
+!
+!*       1.    DEFINITION OF DYF
+!              ------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('DYF',0,ZHOOK_HANDLE)
+!
+CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('DYF',1,ZHOOK_HANDLE)
+END SUBROUTINE DYF_PHY
+!
 END MODULE SHUMAN_PHY
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index f9b556b40bebab63acb1bae728127371e8e14d76..c38f60ff29e27d72d7cc9ecb21f60d0fb09b66b3 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -26,7 +26,7 @@ CONTAINS
 !----------------------------------------------------------------------------
       SUBROUTINE PRANDTL(D,CST,CSTURB,KRR,KSV,KRRI,OTURB_DIAG,&
                          HTURBDIM,OOCEAN,OHARAT,O2D,OCOMPUTE_SRC,&
-                         TPFILE,                               &
+                         TPFILE, OFLAT,                        &
                          PDXX,PDYY,PDZZ,PDZX,PDZY,             &
                          PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
                          PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
@@ -153,6 +153,7 @@ USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_PARAMETERS, ONLY: JPVEXT_TURB
 !
 USE MODI_GRADIENT_M
+USE MODE_GRADIENT_M_PHY, ONLY: GX_M_M_PHY, GY_M_M_PHY
 USE MODE_EMOIST, ONLY : EMOIST
 USE MODE_ETHETA, ONLY : ETHETA
 USE MODI_SHUMAN, ONLY: MZM
@@ -177,6 +178,7 @@ LOGICAL,                INTENT(IN)   ::  OOCEAN       ! switch for Ocean model v
 LOGICAL,                INTENT(IN)   ::  OHARAT
 LOGICAL,                INTENT(IN)   ::  OCOMPUTE_SRC ! flag to define dimensions of SIGS and
 LOGICAL, INTENT(IN) :: O2D               ! Logical for 2D model version (modd_conf)
+LOGICAL,                INTENT(IN)   ::  OFLAT        ! Logical for zero ororography
 CHARACTER(LEN=4),       INTENT(IN)   ::  HTURBDIM     ! Kind of turbulence param.
 TYPE(TFILEDATA),        INTENT(IN)   ::  TPFILE       ! Output file
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)   ::  PDXX,PDYY,PDZZ,PDZX,PDZY
@@ -220,7 +222,8 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT)  ::  PEMOIST ! coefficient E_moi
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT) ::  &
                   ZW1, ZW2, ZW3, &
 !                                                 working variables
-                  ZWORK1,ZWORK2,ZWORK3 
+                  ZWORK1,ZWORK2,ZWORK3,ZWORK4, ZWORK5, ZWORK6,ZWORK7, &
+                  ZGXMM_PTH,ZGYMM_PTH,ZGXMM_PRM,ZGYMM_PRM, ZGXMM_PSV,ZGYMM_PSV
 !                                                 working variables for explicit array
 !                                                     
 INTEGER :: IKB      ! vertical index value for the first inner mass point
@@ -408,81 +411,105 @@ IF(HTURBDIM=='1DIM') THEN        ! 1D case
 !
 ELSE IF (O2D) THEN                      ! 3D case in a 2D model
 !
-    ZWORK1 = MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)**2, D%NKA, D%NKU, D%NKL)
+    CALL GX_M_M_PHY(D,OFLAT,PTHLM,PDXX,PDZZ,PDZX,ZGXMM_PTH)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)**2
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    CALL MZM_PHY(D,ZWORK1,ZWORK2)
+    !
   IF (KRR /= 0) THEN                 ! moist 3D case
-    ZWORK2 = MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)**2, D%NKA, D%NKU, D%NKL)
-    ZWORK3 = MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*     &
-                 GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
-!
-    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2+(CSTURB%XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-                     ZWORK1(:,:,:)
+    CALL GX_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZGXMM_PRM)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)**2
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    CALL MZM_PHY(D,ZWORK1,ZWORK3)
+    !
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    CALL MZM_PHY(D,ZWORK1,ZWORK4)
+    !
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2+(CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                       *PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) )**2 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (CSTURB%XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 * &
-                    ZWORK2(:,:,:)
+    PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2 + (CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                     * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT))**2 * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)
 !
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  CSTURB%XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-                      PEMOIST(:,:,:) * PETHETA(:,:,:) *                         &
-                      ZWORK3(:,:,:)
-    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+    PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) +  CSTURB%XCTV**2 &
+                                        * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)**2   &
+                                        * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                        * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+D%NKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+D%NKL) 
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+D%NKL)
+    PRED2TH3(IIB:IIE,IJB:IJE,IKB)=PRED2TH3(IIB:IIE,IJB:IJE,IKB+D%NKL)
+    PRED2R3(IIB:IIE,IJB:IJE,IKB)=PRED2R3(IIB:IIE,IJB:IJE,IKB+D%NKL) 
+    PRED2THR3(IIB:IIE,IJB:IJE,IKB)=PRED2THR3(IIB:IIE,IJB:IJE,IKB+D%NKL)
 !
   ELSE                 ! dry 3D case in a 2D model
-    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  CSTURB%XCTV**2*PBLL_O_E(:,:,:)**2 *     &
-                      ZWORK1(:,:,:)      
-    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+D%NKL)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) = PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2 +  CSTURB%XCTV**2 & 
+                                       * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)**2 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)      
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2TH3(IIB:IIE,IJB:IJE,IKB)=PRED2TH3(IIB:IIE,IJB:IJE,IKB+D%NKL)
 !
-    PRED2R3(:,:,:) = 0.
+    PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT) = 0.
 !
-    PRED2THR3(:,:,:) = 0.
+    PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT) = 0.
 !
   END IF
 !
 ELSE                                 ! 3D case in a 3D model
 !
-  ZWORK1 = MZM(GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)**2 &
-      + GY_M_M(PTHLM,PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)**2, D%NKA, D%NKU, D%NKL)
-
+  CALL GX_M_M_PHY(D,OFLAT,PTHLM,PDXX,PDZZ,PDZX,ZGXMM_PTH)
+  CALL GY_M_M_PHY(D,OFLAT,PTHLM,PDYY,PDZZ,PDZY,ZGYMM_PTH)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)**2 + ZGYMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)**2
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+  CALL MZM_PHY(D,ZWORK1,ZWORK2)
+  !
   IF (KRR /= 0) THEN                 ! moist 3D case
-    ZWORK2 = MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)**2 + &
-        GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)**2, D%NKA, D%NKU, D%NKL)
-    ZWORK3 = MZM(GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*   &
-         GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)+                           &
-         GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)*                    &
-         GY_M_M(PTHLM,PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL), D%NKA, D%NKU, D%NKL)
-
-    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2TH3(:,:,:)= PREDTH1(:,:,:)**2 +  ( CSTURB%XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:) )**2 * &
-                     ZWORK1(:,:,:)      
+    CALL GX_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZGXMM_PRM)
+    CALL GY_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZGYMM_PRM)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)**2 + ZGYMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)**2
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    CALL MZM_PHY(D,ZWORK1,ZWORK3)
+    !
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    + ZGYMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT) * ZGYMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    CALL MZM_PHY(D,ZWORK1,ZWORK4)
+    !
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2 +  ( CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) )**2 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)      
 !
-    PRED2R3(:,:,:)= PREDR1(:,:,:)**2 + (CSTURB%XCTV*PBLL_O_E(:,:,:)*PEMOIST(:,:,:))**2 *      &
-                    ZWORK2(:,:,:)
+    PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2 + (CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                     * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT))**2 * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)
 !
-    PRED2THR3(:,:,:)= PREDR1(:,:,:) * PREDTH1(:,:,:) +  CSTURB%XCTV**2*PBLL_O_E(:,:,:)**2 *   &
-         PEMOIST(:,:,:) * PETHETA(:,:,:) * ZWORK3(:,:,:)
+    PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)= PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) +  CSTURB%XCTV**2 &
+                                       * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)**2 *   &
+                            PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
 
-    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+D%NKL)
-    PRED2R3(:,:,IKB)=PRED2R3(:,:,IKB+D%NKL)
-    PRED2THR3(:,:,IKB)=PRED2THR3(:,:,IKB+D%NKL)
+    PRED2TH3(IIB:IIE,IJB:IJE,IKB)=PRED2TH3(IIB:IIE,IJB:IJE,IKB+D%NKL)
+    PRED2R3(IIB:IIE,IJB:IJE,IKB)=PRED2R3(IIB:IIE,IJB:IJE,IKB+D%NKL)
+    PRED2THR3(IIB:IIE,IJB:IJE,IKB)=PRED2THR3(IIB:IIE,IJB:IJE,IKB+D%NKL)
 !
   ELSE                 ! dry 3D case in a 3D model
-    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2TH3(:,:,:) = PREDTH1(:,:,:)**2 +  CSTURB%XCTV**2*PBLL_O_E(:,:,:)**2 *                &
-                      ZWORK1(:,:,:)
-    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) = PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2 +  CSTURB%XCTV**2 &
+                                        * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)**2 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ! 
-    PRED2TH3(:,:,IKB)=PRED2TH3(:,:,IKB+D%NKL)
+    PRED2TH3(IIB:IIE,IJB:IJE,IKB)=PRED2TH3(IIB:IIE,IJB:IJE,IKB+D%NKL)
 !
-    PRED2R3(:,:,:) = 0.
+    PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT) = 0.
 !
-    PRED2THR3(:,:,:) = 0.
+    PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT) = 0.
 !
   END IF
 !
@@ -508,76 +535,118 @@ DO JSV=1,KSV
 !
   ELSE  IF (O2D) THEN ! 3D case in a 2D model
 !
-    IF (OOCEAN) THEN    
+    IF (OOCEAN) THEN
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (CST%XG *CST%XALPHAOC * PLM(IIB:IIE,IJB:IJE,1:D%NKT) * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                       / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT))**2
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK2)  
       IF (KRR /= 0) THEN
-        ZW1 = MZM((CST%XG *CST%XALPHAOC * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL) *PETHETA
+        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+        ZW1(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
+        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       ELSE
-        ZW1 = MZM((CST%XG *CST%XALPHAOC * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL)
+        ZW1 = ZWORK2
       END IF
     ELSE
-      ZW1 = MZM((CST%XG / PTHVREF * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL)
-      ZWORK2 = MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*       &
-                             GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL),                 &
-                             D%NKA, D%NKU, D%NKL)
-      ZWORK3 =  MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*       &
-                             GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL),          &
-                             D%NKA, D%NKU, D%NKL)
-!
-      !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * PLM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT))**2
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZW1)  
+      !
+      CALL GX_M_M_PHY(D,OFLAT,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX,ZGXMM_PSV)
+      CALL GX_M_M_PHY(D,OFLAT,PTHLM,PDXX,PDZZ,PDZX,ZGXMM_PTH)
+      CALL GX_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZGXMM_PRM)
+      !
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK2)
+      !
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK3)
+!
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
       IF (KRR /= 0) THEN
-        ZWORK1(:,:,:) = ZW1(:,:,:)*PETHETA(:,:,:)
+        ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZW1(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
       ELSE
-        ZWORK1(:,:,:) = ZW1(:,:,:)
+        ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
       END IF
-      PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                         ZWORK1(:,:,:) * ZWORK2(:,:,:)
+      PRED2THS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)   +        &
+                         ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
                          !
       IF (KRR /= 0) THEN
-        PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                         ZW1(:,:,:) * PEMOIST(:,:,:) * ZWORK3(:,:,:)
+        PRED2RS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)   +        &
+                         ZW1(IIB:IIE,IJB:IJE,1:D%NKT) * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)
       ELSE
-        PRED2RS3(:,:,:,JSV) = 0.
+        PRED2RS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = 0.
       END IF
-      !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     END IF
 !
   ELSE ! 3D case in a 3D model
 !
-    IF (OOCEAN) THEN    
+    IF (OOCEAN) THEN
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (CST%XG *CST%XALPHAOC * PLM(IIB:IIE,IJB:IJE,1:D%NKT) * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                       / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT))**2
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK2)  
       IF (KRR /= 0) THEN
-        ZW1 = MZM((CST%XG *CST%XALPHAOC * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL) *PETHETA
+        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+        ZW1(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
+        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       ELSE
-        ZW1 = MZM((CST%XG *CST%XALPHAOC * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL)
+        ZW1 = ZWORK2
       END IF
-    ELSE 
-      ZW1 = MZM((CST%XG / PTHVREF * PLM * PLEPS / PTKEM)**2, D%NKA, D%NKU, D%NKL)
-      ZWORK2 =  MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*       &
-                             GX_M_M(PTHLM,PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)                  &
-                            +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)*       &
-                             GY_M_M(PTHLM,PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL),                 &
-                             D%NKA, D%NKU, D%NKL)
-      ZWORK3 =  MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)*       &
-                             GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, D%NKA, D%NKU, D%NKL)           &
-                            +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL)*       &
-                             GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY, D%NKA, D%NKU, D%NKL),          &
-                             D%NKA, D%NKU, D%NKL)
-
-      !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+    ELSE
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * PLM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT))**2
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZW1)
+      !
+      CALL GX_M_M_PHY(D,OFLAT,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX,ZGXMM_PSV)
+      CALL GX_M_M_PHY(D,OFLAT,PTHLM,PDXX,PDZZ,PDZX,ZGXMM_PTH)
+      CALL GX_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZGXMM_PRM)
+      CALL GY_M_M_PHY(D,OFLAT,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY,ZGYMM_PSV)
+      CALL GY_M_M_PHY(D,OFLAT,PTHLM,PDYY,PDZZ,PDZY,ZGYMM_PTH)
+      CALL GY_M_M_PHY(D,OFLAT,PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZGYMM_PRM)      
+      !
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      + ZGYMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGYMM_PTH(IIB:IIE,IJB:IJE,1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK2)
+      !
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZGXMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGXMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      + ZGYMM_PSV(IIB:IIE,IJB:IJE,1:D%NKT) * ZGYMM_PRM(IIB:IIE,IJB:IJE,1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      CALL MZM_PHY(D,ZWORK1,ZWORK3)      
+      !
       IF (KRR /= 0) THEN
-        ZWORK1(:,:,:) = ZW1(:,:,:)*PETHETA(:,:,:)
+        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+        ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZW1(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
+        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       ELSE
-        ZWORK1(:,:,:) = ZW1(:,:,:)
+        ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
       END IF
-      PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                         ZWORK1(:,:,:)*ZWORK2(:,:,:)
-                        !
+      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+      PRED2THS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)   +        &
+                         ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       IF (KRR /= 0) THEN
-        PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
-                         ZW1(:,:,:) * PEMOIST(:,:,:) * ZWORK3(:,:,:)
+        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+        PRED2RS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)   +        &
+                         ZW1(IIB:IIE,IJB:IJE,1:D%NKT) * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)
+        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       ELSE
-        PRED2RS3(:,:,:,JSV) = 0.
+        PRED2RS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = 0.
       END IF
-      !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
   END IF 
 !
   END IF ! end of HTURBDIM if-block
@@ -2127,7 +2196,7 @@ IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)&
 * (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2          &
-        +(1.+2.*PREDR1(:,:,:))/PD(:,:,:)
+        +(1.+2.*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
diff --git a/src/common/turb/mode_turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
index 2df64d99780f12d7e692fdfc9f178d754bc0fc22..49044ab5c2f39f86710b37b43023f2aea47b9f0f 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -419,7 +419,7 @@ IJB=D%NJBC
 !
 CALL PRANDTL(D,CST,CSTURB,KRR,KSV,KRRI,OTURB_FLX,  &
              HTURBDIM,OOCEAN,OHARAT,O2D,OCOMPUTE_SRC,&
-             TPFILE,                               &
+             TPFILE, OFLAT,                        &
              PDXX,PDYY,PDZZ,PDZX,PDZY,             &
              PTHVREF,PLOCPEXNM,PATHETA,PAMOIST,    &
              PLM,PLEPS,PTKEM,PTHLM,PRM,PSVM,PSRCM, &
diff --git a/src/mesonh/aux/shuman_phy.f90 b/src/mesonh/aux/shuman_phy.f90
index 368ea93330417abb96e4e77fe6aa238e4cd41938..fa3f78966eeaf940d22d24a7b41ba8947e67912b 100644
--- a/src/mesonh/aux/shuman_phy.f90
+++ b/src/mesonh/aux/shuman_phy.f90
@@ -786,4 +786,213 @@ END DO
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE DZF_PHY
+!
+!     ###############################
+      SUBROUTINE DXF_PHY(D,PA,PDXF)
+!     ###############################
+!
+!!****  *DXF* -  Shuman operator : finite difference operator in x direction
+!!                                  for a variable at a flux side
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function  is to compute a finite difference 
+!     along the x direction (I index) for a field PA localized at a x-flux
+!     point (u point). The result is localized at a mass point.
+!
+!!**  METHOD
+!!    ------ 
+!!        The result PDXF(i,:,:) is defined by (PA(i+1,:,:)-PA(i,:,:))
+!!        At i=size(PA,1), PDXF(i,:,:) are replaced by the values of PDXF,
+!!    which are the right values in the x-cyclic case
+!!    
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_PARAMETERS: declaration of parameter variables
+!!        JPHEXT: define the number of marginal points out of the 
+!!        physical domain along the horizontal directions.
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (SHUMAN operators)
+!!      Technical specifications Report of The Meso-NH (chapters 3)  
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!	V. Ducrocq       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/07/94 
+!!      Modification to include the periodic case 13/10/94 J.Stein 
+!!                   optimisation                 20/08/00 J. Escobar
+!!      correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1    
+!!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS, ONLY: JPHEXT
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of argument and result
+!              ------------------------------------
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)                :: PA     ! variable at flux
+                                                            !  side
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PDXF   ! result at mass
+                                                            ! localization 
+!
+!*       0.2   Declarations of local variables
+!              -------------------------------
+!
+INTEGER :: JI             ! Loop index in x direction
+INTEGER :: IIU            ! upper bound in x direction of PA 
+!             
+INTEGER :: JJK,IJU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!             
+!
+!-------------------------------------------------------------------------------
+!
+!*       1.    DEFINITION OF DXF
+!              ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR  = 1 + 1
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+   PDXF(JIJK-1,1,1) = PA(JIJK,1,1) - PA(JIJK-1,1,1) 
+END DO
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JI=1,JPHEXT
+   DO JJK=1,IJU*IKU
+      PDXF(IIU-JPHEXT+JI,JJK,1) = PDXF(JPHEXT+JI,JJK,1) ! for reprod JPHEXT <> 1
+   END DO
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE DXF_PHY
+!
+!     ###############################
+      SUBROUTINE DYF_PHY(D,PA,PDYF)
+!     ###############################
+!
+!!****  *DYF* -  Shuman operator : finite difference operator in y direction
+!!                                  for a variable at a flux side
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function  is to compute a finite difference 
+!     along the y direction (J index) for a field PA localized at a y-flux
+!     point (v point). The result is localized at a mass point.
+!
+!!**  METHOD
+!!    ------ 
+!!        The result PDYF(:,j,:) is defined by (PA(:,j+1,:)-PA(:,j,:))
+!!        At j=size(PA,2), PDYF(:,j,:) are replaced by the values of PDYM,
+!!    which are the right values in the y-cyclic case
+!!    
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_PARAMETERS: declaration of parameter variables
+!!        JPHEXT: define the number of marginal points out of the 
+!!        physical domain along the horizontal directions.
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH (SHUMAN operators)
+!!      Technical specifications Report of The Meso-NH (chapters 3)  
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!	V. Ducrocq       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/07/94 
+!!      Modification to include the periodic case 13/10/94 J.Stein 
+!!                   optimisation                 20/08/00 J. Escobar
+!!      correction of in halo/pseudo-cyclic calculation for JPHEXT<> 1 
+!!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_PARAMETERS, ONLY: JPHEXT
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of argument and result
+!              ------------------------------------
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)  :: D
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)                :: PA     ! variable at flux
+                                                            !  side
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PDYF   ! result at mass
+                                                            ! localization 
+!
+!*       0.2   Declarations of local variables
+!              -------------------------------
+!
+INTEGER :: JJ            ! Loop index in y direction
+INTEGER :: IJU           ! upper bound in y direction of PA 
+!
+!          
+INTEGER :: IIU,IKU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+!            
+!-------------------------------------------------------------------------------
+!
+!*       1.    DEFINITION OF DYF
+!              ------------------
+!
+IIU = SIZE(PA,1)
+IJU = SIZE(PA,2)
+IKU = SIZE(PA,3)
+!
+JIJKOR  = 1 + IIU
+JIJKEND = IIU*IJU*IKU
+!
+!CDIR NODEP
+!OCL NOVREC
+DO JIJK=JIJKOR , JIJKEND
+   PDYF(JIJK-IIU,1,1)         = PA(JIJK,1,1)  -  PA(JIJK-IIU,1,1) 
+END DO
+!
+DO JJ=1,JPHEXT
+   PDYF(:,IJU-JPHEXT+JJ,:) = PDYF(:,JPHEXT+JJ,:) ! for reprod JPHEXT <> 1
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE DYF_PHY
+!
 END MODULE SHUMAN_PHY