From b02565901b9e6485abef2e2dda2a2061cfcdba4e Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 3 May 2022 12:54:46 +0200
Subject: [PATCH] Quentin 03/05/2022: Expand in physical points : mode_prandtl
 (only turb 1D) + gradient GZ_M_W at physical points and in subroutine

---
 src/common/aux/gradient_m_phy.F90 |  102 +++
 src/common/turb/mode_prandtl.F90  | 1251 +++++++++++++++++------------
 2 files changed, 822 insertions(+), 531 deletions(-)
 create mode 100644 src/common/aux/gradient_m_phy.F90

diff --git a/src/common/aux/gradient_m_phy.F90 b/src/common/aux/gradient_m_phy.F90
new file mode 100644
index 000000000..fb2119e05
--- /dev/null
+++ b/src/common/aux/gradient_m_phy.F90
@@ -0,0 +1,102 @@
+MODULE MODE_GRADIENT_M_PHY
+IMPLICIT NONE
+CONTAINS
+!     #########################################
+      SUBROUTINE GZ_M_W_PHY(D,PY,PDZZ,PGZ_M_W)
+!     #########################################
+!
+!!****  *GZ_M_W * - Compute the gradient along z direction for a 
+!!       variable localized at a mass point
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this routine is to compute a gradient along x,y,z 
+!     directions for a field PY localized at a mass point. The result PGZ_M_W
+!     is localized at a z-flux point (w point)
+!
+!              
+!                    dzm(PY)  
+!       PGZ_M_W =    ------- 
+!                     d*zz        
+!
+!!**  METHOD
+!!    ------
+!!      We employ the Shuman operators to compute the derivatives and the 
+!!    averages. The metric coefficients PDZZ are dummy arguments.
+!!
+!!
+!!    EXTERNAL
+!!    --------
+!!      FUNCTION DZM : compute a finite difference along the z 
+!!    direction for a variable at a mass localization
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------  
+!!      Module MODI_SHUMAN : interface for the Shuman functions
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation (function GZ_M_W)
+!!      
+!!
+!!    AUTHOR
+!!    ------
+!!	P. Hereil and J. Stein       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    05/07/94 
+!!      Modification       16/03/95  change the order of the arguments
+!!                         19/07/00  inlining(J. Stein)
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of arguments and results
+!              -------------------------------------
+!
+TYPE(DIMPHYEX_t),       INTENT(IN)   :: D
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)  :: PDZZ                   !d*zz
+!
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)                :: PY       ! variable at mass
+                                                              ! localization
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(OUT) :: PGZ_M_W  ! result at flux
+                                                              ! side
+!
+INTEGER :: IKT,IKTB,IKTE,IIB,IJB,IIE,IJE
+INTEGER :: JI,JJ,JK
+!-------------------------------------------------------------------------------
+!
+!*       1.    COMPUTE THE GRADIENT ALONG Z
+!              -----------------------------
+!
+IKT=D%NKT
+IKTB=D%NKTB              
+IKTE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+DO JK=IKTB,IKTE 
+  DO JJ=IJB,IJE 
+    DO JI=IIB,IIE 
+      PGZ_M_W(JI,JJ,JK) =  (PY(JI,JJ,JK)-PY(JI,JJ,JK-D%NKL )) / PDZZ(JI,JJ,JK)
+  ENDDO
+ ENDDO
+ENDDO
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+PGZ_M_W(IIB:IIE,IJB:IJE,D%NKU)=  (PY(IIB:IIE,IJB:IJE,D%NKU)-PY(IIB:IIE,IJB:IJE,D%NKU-D%NKL))  &
+                           / PDZZ(IIB:IIE,IJB:IJE,D%NKU)
+PGZ_M_W(IIB:IIE,IJB:IJE,D%NKA)= PGZ_M_W(IIB:IIE,IJB:IJE,D%NKU) ! -999.
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE GZ_M_W_PHY
+END MODULE MODE_GRADIENT_M_PHY
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index 2ba414d71..6c391103f 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -18,6 +18,8 @@ USE MODD_DIMPHYEX,   ONLY : DIMPHYEX_t
 USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
 !
 USE MODI_SHUMAN, ONLY: MZM, MZF
+USE SHUMAN_PHY, ONLY: MZM_PHY,MZF_PHY
+USE MODE_GRADIENT_M_PHY
 IMPLICIT NONE
 !----------------------------------------------------------------------------
 CONTAINS
@@ -223,7 +225,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) ::  &
 !                                                     
 INTEGER :: IKB      ! vertical index value for the first inner mass point
 INTEGER :: IKE      ! vertical index value for the last inner mass point
-INTEGER::  JSV,JI,JJ,JK ! loop index
+INTEGER::  JSV,JI,JJ,JK,IIB,IIE,IJB,IJE ! loop index
 
 INTEGER :: JLOOP
 REAL    :: ZMINVAL
@@ -248,34 +250,44 @@ PRED2RS3(:,:,:,:)=0.
 PBLL_O_E(:,:,:)=0.
 ENDIF
 !
-IKB = D%NKB
-IKE = D%NKE 
-!
-PETHETA(:,:,:) = MZM(ETHETA(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM,OOCEAN,OCOMPUTE_SRC), D%NKA, D%NKU, D%NKL)
-PEMOIST(:,:,:) = MZM(EMOIST(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM,OOCEAN), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
-PETHETA(:,:,D%NKA) = 2.*PETHETA(:,:,IKB) - PETHETA(:,:,IKB+D%NKL)
-PEMOIST(:,:,D%NKA) = 2.*PEMOIST(:,:,IKB) - PEMOIST(:,:,IKB+D%NKL)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+IKB=D%NKTB
+IKE=D%NKTE 
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+ZWORK1 = ETHETA(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM,OOCEAN,OCOMPUTE_SRC)
+ZWORK2 = EMOIST(D,CST,KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM,OOCEAN)
+CALL MZM_PHY(D,ZWORK1,PETHETA)
+CALL MZM_PHY(D,ZWORK2,PEMOIST)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
+PETHETA(IIB:IIE,IJB:IJE,D%NKA) = 2.*PETHETA(IIB:IIE,IJB:IJE,IKB) - PETHETA(IIB:IIE,IJB:IJE,IKB+D%NKL)
+PEMOIST(IIB:IIE,IJB:IJE,D%NKA) = 2.*PEMOIST(IIB:IIE,IJB:IJE,IKB) - PEMOIST(IIB:IIE,IJB:IJE,IKB+D%NKL)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 !
 !---------------------------------------------------------------------------
 IF (.NOT. OHARAT) THEN
 !
 !          1.3 1D Redelsperger numbers
 !
-PBLL_O_E(:,:,:) = MZM(CST%XG / PTHVREF(:,:,:) * PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:), D%NKA, D%NKU, D%NKL)
-ZWORK1 =  GZ_M_W(D%NKA, D%NKU, D%NKL,PTHLM,PDZZ)
+!$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)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+CALL MZM_PHY(D,ZWORK1,PBLL_O_E)
+!
+CALL GZ_M_W_PHY(D,PTHLM,PDZZ,ZWORK1)
 IF (KRR /= 0) THEN                ! moist case
-  ZWORK2 =  GZ_M_W(D%NKA, D%NKU, D%NKL,PRM(:,:,:,1),PDZZ)
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  PREDTH1(:,:,:)= CSTURB%XCTV*PBLL_O_E(:,:,:) * PETHETA(:,:,:) * ZWORK1(:,:,:)
-  PREDR1(:,:,:) = CSTURB%XCTV*PBLL_O_E(:,:,:) * PEMOIST(:,:,:) * ZWORK2(:,:,:)
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+  CALL GZ_M_W_PHY(D,PRM(:,:,:,1),PDZZ,ZWORK2)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)= CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+  PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PEMOIST(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)
 ELSE                              ! dry case
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  PREDTH1(:,:,:)= CSTURB%XCTV*PBLL_O_E(:,:,:)  * ZWORK1(:,:,:)
-  PREDR1(:,:,:) = 0.
-  !$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)
+  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)= CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)  * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+  PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) = 0.
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 END IF
 !
 !       3. Limits on 1D Redelperger numbers
@@ -284,8 +296,8 @@ END IF
 ZMINVAL = (1.-1./CSTURB%XPHI_LIM)
 !
 DO JK=1,D%NKT 
- DO JJ=1,D%NJT 
-  DO JI=1,D%NIT 
+ DO JJ=IJB,IJE
+  DO JI=IIB,IIE 
    ZW1(JI,JJ,JK) = 1.
    ZW2(JI,JJ,JK) = 1.
    !
@@ -339,16 +351,16 @@ ENDDO
 !
 !          For the scalar variables
 DO JSV=1,KSV
-  ZWORK1 = GZ_M_W(D%NKA, D%NKU, D%NKL,PSVM(:,:,:,JSV),PDZZ)
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  PREDS1(:,:,:,JSV)=CSTURB%XCTV*PBLL_O_E(:,:,:)*ZWORK1(:,:,:)
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+  CALL GZ_M_W_PHY(D,PSVM(:,:,:,JSV),PDZZ,ZWORK1)
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+  PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)=CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 END DO
 !
 DO JSV=1,KSV
  DO JK=1,D%NKT 
-  DO JJ=1,D%NJT 
-   DO JI=1,D%NIT 
+  DO JJ=IJB,IJE
+   DO JI=IIB,IIE 
     IF(PREDS1(JI,JJ,JK,JSV) < 0.) THEN
      ZW2(JI,JJ,JK)=-1.
     ELSE
@@ -368,13 +380,13 @@ ENDDO
 IF(HTURBDIM=='1DIM') THEN        ! 1D case
 !
 !
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  PRED2TH3(:,:,:)  = PREDTH1(:,:,:)**2
+  !$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
 !
-  PRED2R3(:,:,:)   = PREDR1(:,:,:) **2
+  PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT)   = PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) **2
 !
-  PRED2THR3(:,:,:) = PREDTH1(:,:,:) * PREDR1(:,:,:)
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+  PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT) = PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 ELSE IF (O2D) THEN                      ! 3D case in a 2D model
 !
@@ -467,14 +479,14 @@ DO JSV=1,KSV
 !
   IF(HTURBDIM=='1DIM') THEN
 !        1D case
-    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-    PRED2THS3(:,:,:,JSV)  = PREDS1(:,:,:,JSV) * PREDTH1(:,:,:)
+    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+    PRED2THS3(IIB:IIE,IJB:IJE,1:D%NKT,JSV)  = PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV) * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)
     IF (KRR /= 0) THEN
-      PRED2RS3(:,:,:,JSV)   = PREDR1(:,:,:) *PREDS1(:,:,:,JSV)
+      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)
     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)
 !
   ELSE  IF (O2D) THEN ! 3D case in a 2D model
 !
@@ -628,20 +640,24 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN)    :: PF_LIM  ! Value of F when P
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PF      ! function F to smooth
 !
 REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZCOEF
-INTEGER :: JI,JJ,JK
+INTEGER :: JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 !* adds a artificial correction to smooth the function near the discontinuity
 !  point at Phi3 = Phi_lim
 !  This smoothing is applied between 0.9*phi_lim (=2.7) and Phi_lim (=3)
 !   Note that in the Boundary layer, phi is usually between 0.8 and 1
 !
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 !
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
-ZCOEF(:,:,:) = MAX(MIN((  10.*(1.-PPHI3(:,:,:)/CSTURB%XPHI_LIM)) ,1.), 0.) 
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+ZCOEF(IIB:IIE,IJB:IJE,1:D%NKT) = MAX(MIN((  10.*(1.-PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XPHI_LIM)) ,1.), 0.) 
 !
-PF(:,:,:) =     ZCOEF(:,:,:)   * PF(:,:,:)    &
-          + (1.-ZCOEF(:,:,:))  * PF_LIM(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+PF(IIB:IIE,IJB:IJE,1:D%NKT) =     ZCOEF(IIB:IIE,IJB:IJE,1:D%NKT)   * PF(IIB:IIE,IJB:IJE,1:D%NKT)    &
+          + (1.-ZCOEF(IIB:IIE,IJB:IJE,1:D%NKT))  * PF_LIM(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 !
 END SUBROUTINE SMOOTH_TURB_FUNCT
 !----------------------------------------------------------------------------
@@ -658,58 +674,62 @@ FUNCTION PHI3(D,CSTURB,PREDTH1,PREDR1,PRED2TH3,PRED2R3,PRED2THR3,HTURBDIM,OUSERV
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: PHI3
 !
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZW1, ZW2
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:PHI3',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 !
 IF (HTURBDIM=='3DIM') THEN
         !* 3DIM case
-  !$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)    
   IF (OUSERV) THEN
-    ZW1(:,:,:) = 1. + 1.5* (PREDTH1(:,:,:)+PREDR1(:,:,:)) +      &
-                   ( 0.5 * (PREDTH1(:,:,:)**2+PREDR1(:,:,:)**2)  &
-                         + PREDTH1(:,:,:) * PREDR1(:,:,:)        &
+    ZW1(IIB:IIE,IJB:IJE,1:D%NKT) = 1. + 1.5* (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) +      &
+                   ( 0.5 * (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)  &
+                         + PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)        &
                    )
 
-    ZW2(:,:,:) = 0.5 * (PRED2TH3(:,:,:)-PRED2R3(:,:,:))
+    ZW2(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5 * (PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)-PRED2R3(IIB:IIE,IJB:IJE,1:D%NKT))
 
-    PHI3(:,:,:)= 1. -                                          &
-    ( ( (1.+PREDR1(:,:,:)) *                                   &
-        (PRED2THR3(:,:,:) + PRED2TH3(:,:,:)) / PREDTH1(:,:,:)  &
-      ) + ZW2(:,:,:)                                           &
-    ) / ZW1(:,:,:)
+    PHI3(IIB:IIE,IJB:IJE,1:D%NKT)= 1. -                                          &
+    ( ( (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) *                                   &
+        (PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT) + PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)) / PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)  &
+      ) + ZW2(IIB:IIE,IJB:IJE,1:D%NKT)                                           &
+    ) / ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
   ELSE
-    ZW1(:,:,:) = 1. + 1.5* PREDTH1(:,:,:) + &
-                 0.5* PREDTH1(:,:,:)**2
+    ZW1(IIB:IIE,IJB:IJE,1:D%NKT) = 1. + 1.5* PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + &
+                 0.5* PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)**2
 
-    ZW2(:,:,:) = 0.5* PRED2TH3(:,:,:)
+    ZW2(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5* PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)
 
-    PHI3(:,:,:)= 1. -                                       &
-            (PRED2TH3(:,:,:) / PREDTH1(:,:,:) + ZW2(:,:,:)) / ZW1(:,:,:)
+    PHI3(IIB:IIE,IJB:IJE,1:D%NKT)= 1. -                                       &
+            (PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) / PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + ZW2(IIB:IIE,IJB:IJE,1:D%NKT)) / ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
   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)    
 
-  !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
-  WHERE( PHI3(:,:,:) <= 0. .OR. PHI3(:,:,:) > CSTURB%XPHI_LIM )
-    PHI3(:,:,:) = CSTURB%XPHI_LIM
+  !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  WHERE( PHI3(IIB:IIE,IJB:IJE,1:D%NKT) <= 0. .OR. PHI3(IIB:IIE,IJB:IJE,1:D%NKT) > CSTURB%XPHI_LIM )
+    PHI3(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XPHI_LIM
   END WHERE
-  !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+  !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 ELSE
         !* 1DIM case
-  !$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)
   IF (OUSERV) THEN
-    PHI3(:,:,:)= 1./(1.+PREDTH1(:,:,:)+PREDR1(:,:,:))
+    PHI3(IIB:IIE,IJB:IJE,1:D%NKT)= 1./(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))
   ELSE
-    PHI3(:,:,:)= 1./(1.+PREDTH1(:,:,:))
+    PHI3(IIB:IIE,IJB:IJE,1:D%NKT)= 1./(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))
   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
 !
-PHI3(:,:,IKB-1)=PHI3(:,:,IKB)
-PHI3(:,:,IKE+1)=PHI3(:,:,IKE)
+PHI3(IIB:IIE,IJB:IJE,IKB-1)=PHI3(IIB:IIE,IJB:IJE,IKB)
+PHI3(IIB:IIE,IJB:IJE,IKE+1)=PHI3(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:PHI3',1,ZHOOK_HANDLE)
 END FUNCTION PHI3
@@ -727,36 +747,40 @@ FUNCTION PSI_SV(D,CSTURB,KSV,PREDTH1,PREDR1,PREDS1,PRED2THS,PRED2RS,PPHI3,PPSI3)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT),   INTENT(IN) :: PPSI3
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KSV) :: PSI_SV
 !
-  INTEGER :: IKB, IKE
+  INTEGER :: IKB, IKE, IIB,IIE,IJB,IJE
   INTEGER :: JSV,JI,JJ,JK
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:PSI_SV',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
-!
-DO JSV=1,SIZE(PSI_SV,4)
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
-  PSI_SV(:,:,:,JSV) = ( 1.                                             &
-    - (CSTURB%XCPR3+CSTURB%XCPR5) * (PRED2THS(:,:,:,JSV)/PREDS1(:,:,:,JSV)-PREDTH1(:,:,:)) &
-    - (CSTURB%XCPR4+CSTURB%XCPR5) * (PRED2RS(:,:,:,JSV)/PREDS1(:,:,:,JSV)-PREDR1(:,:,:)) &
-    - CSTURB%XCPR3 * PREDTH1(:,:,:) * PPHI3(:,:,:) - CSTURB%XCPR4 * PREDR1(:,:,:) * PPSI3(:,:,:)                 &
-                ) / (  1. + CSTURB%XCPR5 * ( PREDTH1(:,:,:) + PREDR1(:,:,:) ) )           
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!
+DO JSV=1,KSV
+  !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = ( 1.                                             &
+    - (CSTURB%XCPR3+CSTURB%XCPR5) * (PRED2THS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+    - (CSTURB%XCPR4+CSTURB%XCPR5) * (PRED2RS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+    - CSTURB%XCPR3 * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) - CSTURB%XCPR4 * PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PPSI3(IIB:IIE,IJB:IJE,1:D%NKT)                 &
+                ) / (  1. + CSTURB%XCPR5 * ( PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + PREDR1(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)    
 !        control of the PSI_SV positivity
-  !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
-  WHERE ( (PSI_SV(:,:,:,JSV) <=0.).AND. (PREDTH1(:,:,:)+PREDR1(:,:,:)) <= 0. )
-    PSI_SV(:,:,:,JSV)=CSTURB%XPHI_LIM
+  !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
+  WHERE ( (PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) <=0.).AND. (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) <= 0. )
+    PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV)=CSTURB%XPHI_LIM
   END WHERE
-  !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+  !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
-  PSI_SV(:,:,:,JSV) = MAX( 1.E-4, MIN(CSTURB%XPHI_LIM,PSI_SV(:,:,:,JSV)) )
-  !$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)    
+  PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = MAX( 1.E-4, MIN(CSTURB%XPHI_LIM,PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV)) )
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 !
-  PSI_SV(:,:,IKB-1,JSV)=PSI_SV(:,:,IKB,JSV)
-  PSI_SV(:,:,IKE+1,JSV)=PSI_SV(:,:,IKE,JSV)
+  PSI_SV(IIB:IIE,IJB:IJE,IKB-1,JSV)=PSI_SV(IIB:IIE,IJB:IJE,IKB,JSV)
+  PSI_SV(IIB:IIE,IJB:IJE,IKE+1,JSV)=PSI_SV(IIB:IIE,IJB:IJE,IKE,JSV)
 END DO
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:PSI_SV',1,ZHOOK_HANDLE)
@@ -773,64 +797,68 @@ FUNCTION D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTU
   CHARACTER(LEN=4),       INTENT(IN) :: HTURBDIM  ! 1DIM or 3DIM turb. scheme
   LOGICAL,                INTENT(IN) :: OUSERV    ! flag to use vapor
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_PHI3DTDZ_O_DDTDZ
-  INTEGER :: IKB, IKE,JL,JK,JJ,JI
+  INTEGER :: IKB, IKE,JK,JJ,JI, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DTDZ_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 !
 IF (HTURBDIM=='3DIM') THEN
         !* 3DIM case
   IF (OUSERV) THEN
-   !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 #ifdef REPRO48
-    WHERE (PPHI3(:,:,:)/=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/=CSTURB%XPHI_LIM)
 #else
-    WHERE (PPHI3(:,:,:)<=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-      D_PHI3DTDZ_O_DDTDZ(:,:,:) = PPHI3(:,:,:)                       &
-          * (1. - PREDTH1(:,:,:) * (3./2.+PREDTH1(:,:,:)+PREDR1(:,:,:))          &
-               /((1.+PREDTH1(:,:,:)+PREDR1(:,:,:))*(1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:))))) &
-          + (1.+PREDR1(:,:,:))*(PRED2THR3(:,:,:)+PRED2TH3(:,:,:))                       &
-               / (PREDTH1(:,:,:)*(1.+PREDTH1(:,:,:)+PREDR1(:,:,:))* &
-                 (1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:)))) &
-          - (1./2.*PREDTH1(:,:,:)+PREDR1(:,:,:) * (1.+PREDTH1(:,:,:)+PREDR1(:,:,:)))           &
-               / ((1.+PREDTH1(:,:,:)+PREDR1(:,:,:))*(1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:))))
+      D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)                       &
+          * (1. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * (3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))          &
+               /((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))) &
+          + (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT))                       &
+               / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))* &
+                 (1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))) &
+          - (1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))           &
+               / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
     ELSEWHERE
-      D_PHI3DTDZ_O_DDTDZ(:,:,:) = PPHI3(:,:,:)
+      D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     ENDWHERE
-   !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 
 !
   ELSE
-   !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 #ifdef REPRO48
-    WHERE (PPHI3(:,:,:)/=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/=CSTURB%XPHI_LIM)
 #else
-    WHERE (PPHI3(:,:,:)<=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-    D_PHI3DTDZ_O_DDTDZ(:,:,:) = PPHI3(:,:,:)             &
-          * (1. - PREDTH1(:,:,:) * (3./2.+PREDTH1(:,:,:))      &
-               /((1.+PREDTH1(:,:,:))*(1.+1./2.*PREDTH1(:,:,:))))        &
-          + PRED2TH3(:,:,:) / (PREDTH1(:,:,:)*(1.+PREDTH1(:,:,:))*(1.+1./2.*PREDTH1(:,:,:))) &
-          - 1./2.*PREDTH1(:,:,:) / ((1.+PREDTH1(:,:,:))*(1.+1./2.*PREDTH1(:,:,:)))
+    D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)             &
+          * (1. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * (3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))      &
+               /((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))))        &
+          + PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))) &
+          - 1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)))
     ELSEWHERE
-      D_PHI3DTDZ_O_DDTDZ(:,:,:) = PPHI3(:,:,:)
+      D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     ENDWHERE
-   !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 !
   END IF
 ELSE
         !* 1DIM case
-DO JJ=1,SIZE(PPHI3,2)
-  DO JL=1,SIZE(PPHI3,1)
-    DO JK=1,SIZE(PPHI3,3)
-      IF ( ABS(PPHI3(JL,JJ,JK)-CSTURB%XPHI_LIM) < 1.E-12 ) THEN
-         D_PHI3DTDZ_O_DDTDZ(JL,JJ,JK)=PPHI3(JL,JJ,JK)*&
-&       (1. - PREDTH1(JL,JJ,JK)*PPHI3(JL,JJ,JK))
+DO JK=1,D%NKT
+  DO JJ=IJB,IJE
+    DO JI=IIB,IIE
+      IF ( ABS(PPHI3(JI,JJ,JK)-CSTURB%XPHI_LIM) < 1.E-12 ) THEN
+         D_PHI3DTDZ_O_DDTDZ(JI,JJ,JK)=PPHI3(JI,JJ,JK)*&
+&       (1. - PREDTH1(JI,JJ,JK)*PPHI3(JI,JJ,JK))
       ELSE
-         D_PHI3DTDZ_O_DDTDZ(JL,JJ,JK)=PPHI3(JL,JJ,JK)
+         D_PHI3DTDZ_O_DDTDZ(JI,JJ,JK)=PPHI3(JI,JJ,JK)
       ENDIF
     ENDDO
   ENDDO
@@ -843,8 +871,8 @@ END IF
 CALL SMOOTH_TURB_FUNCT(D,CSTURB,PPHI3,PPHI3,D_PHI3DTDZ_O_DDTDZ)
 #endif
 !
-D_PHI3DTDZ_O_DDTDZ(:,:,IKB-1)=D_PHI3DTDZ_O_DDTDZ(:,:,IKB)
-D_PHI3DTDZ_O_DDTDZ(:,:,IKE+1)=D_PHI3DTDZ_O_DDTDZ(:,:,IKE)
+D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DTDZ_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_PHI3DTDZ_O_DDTDZ
@@ -860,51 +888,55 @@ FUNCTION D_PHI3DRDZ_O_DDRDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTU
   CHARACTER(LEN=4),       INTENT(IN) :: HTURBDIM  ! 1DIM or 3DIM turb. scheme
   LOGICAL,                INTENT(IN) :: OUSERV    ! flag to use vapor
    REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_PHI3DRDZ_O_DDRDZ
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DRDZ_O_DDRDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 !
 !
 IF (HTURBDIM=='3DIM') THEN
         !* 3DIM case
   IF (OUSERV) THEN
-   !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 #ifdef REPRO48
-    WHERE (PPHI3(:,:,:)/=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/=CSTURB%XPHI_LIM)
 #else
-    WHERE (PPHI3(:,:,:)<=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-      D_PHI3DRDZ_O_DDRDZ(:,:,:) =          &
-                PPHI3(:,:,:) * (1.-PREDR1(:,:,:)*(3./2.+PREDTH1(:,:,:)+PREDR1(:,:,:)) &
-                  / ((1.+PREDTH1(:,:,:)+PREDR1(:,:,:))*(1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:)))))  &
-              - PREDR1(:,:,:) * (PRED2THR3(:,:,:)+PRED2TH3(:,:,:)) / (PREDTH1(:,:,:)         &
-                  * (1.+PREDTH1(:,:,:)+PREDR1(:,:,:))*(1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:))))    &
-              + PREDR1(:,:,:) * (1./2.+PREDTH1(:,:,:)+PREDR1(:,:,:))                  &
-                  / ((1.+PREDTH1(:,:,:)+PREDR1(:,:,:))*(1.+1./2.*(PREDTH1(:,:,:)+PREDR1(:,:,:))))
+      D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) =          &
+                PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) * (1.-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+                  / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))))  &
+              - PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)) / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)         &
+                  * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))    &
+              + PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (1./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))                  &
+                  / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
     ELSEWHERE
-      D_PHI3DRDZ_O_DDRDZ(:,:,:) = PPHI3(:,:,:)
+      D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     END WHERE
-   !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+   !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
   ELSE
-    D_PHI3DRDZ_O_DDRDZ(:,:,:) = PPHI3(:,:,:)
+    D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
   END IF
 ELSE
         !* 1DIM case
-    !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+    !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 #ifdef REPRO48
-    WHERE (PPHI3(:,:,:)/=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/=CSTURB%XPHI_LIM)
 #else
-    WHERE (PPHI3(:,:,:)<=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-    D_PHI3DRDZ_O_DDRDZ(:,:,:) = PPHI3(:,:,:)                           &
-          * (1. - PREDR1(:,:,:)*PPHI3(:,:,:))
+    D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)                           &
+          * (1. - PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*PPHI3(IIB:IIE,IJB:IJE,1:D%NKT))
   ELSEWHERE
-    D_PHI3DRDZ_O_DDRDZ(:,:,:) = PPHI3(:,:,:)
+    D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
   END WHERE
-  !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+  !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 END IF
 !
 #ifdef REPRO48
@@ -913,8 +945,8 @@ END IF
 CALL SMOOTH_TURB_FUNCT(D,CSTURB,PPHI3,PPHI3,D_PHI3DRDZ_O_DDRDZ)
 #endif
 !
-D_PHI3DRDZ_O_DDRDZ(:,:,IKB-1)=D_PHI3DRDZ_O_DDRDZ(:,:,IKB)
-D_PHI3DRDZ_O_DDRDZ(:,:,IKE+1)=D_PHI3DRDZ_O_DDRDZ(:,:,IKE)
+D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,IKB-1)=D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,IKB)
+D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,IKE+1)=D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DRDZ_O_DDRDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_PHI3DRDZ_O_DDRDZ
@@ -932,34 +964,38 @@ FUNCTION D_PHI3DTDZ2_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,PD
   LOGICAL,                INTENT(IN) :: OUSERV    ! flag to use vapor
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_PHI3DTDZ2_O_DDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DTDZ2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 !
 !
 IF (HTURBDIM=='3DIM') THEN
    ! by derivation of (phi3 dtdz) * dtdz according to dtdz we obtain:
    ZWORK1 = D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,OUSERV) 
-   !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-   D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PDTDZ(:,:,:) * (PPHI3(:,:,:) +  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)
+   D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) +  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT))
+   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
         !* 1DIM case
-    !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+    !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 #ifdef REPRO48
-    WHERE (PPHI3(:,:,:)/=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)/=CSTURB%XPHI_LIM)
 #else
-    WHERE (PPHI3(:,:,:)<=CSTURB%XPHI_LIM)
+    WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:)*PDTDZ(:,:,:)             &
-          * (2. - PREDTH1(:,:,:)*PPHI3(:,:,:))
+      D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)             &
+          * (2. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*PPHI3(IIB:IIE,IJB:IJE,1:D%NKT))
     ELSEWHERE
-      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:) * 2. * PDTDZ(:,:,:)
+      D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) * 2. * PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)
     END WHERE
-    !$mnh_end_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
+    !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 END IF
 !
 #ifdef REPRO48
@@ -969,8 +1005,8 @@ CALL SMOOTH_TURB_FUNCT(D,CSTURB,PPHI3,PPHI3*2.*PDTDZ,D_PHI3DTDZ2_O_DDTDZ)
 #endif
 !
 !
-D_PHI3DTDZ2_O_DDTDZ(:,:,IKB-1)=D_PHI3DTDZ2_O_DDTDZ(:,:,IKB)
-D_PHI3DTDZ2_O_DDTDZ(:,:,IKE+1)=D_PHI3DTDZ2_O_DDTDZ(:,:,IKE)
+D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_PHI3DTDZ2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_PHI3DTDZ2_O_DDTDZ
@@ -984,19 +1020,23 @@ FUNCTION M3_WTH_WTH2(D,CSTURB,PREDTH1,PREDR1,PD,PBLL_O_E,PETHETA)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_WTH_WTH2
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WTH2',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_WTH2(:,:,:) = CSTURB%XCSHF*PBLL_O_E(:,:,:)*PETHETA(:,:,:)*0.5/CSTURB%XCTD        &
-                   * (1.+0.5*PREDTH1(:,:,:)+PREDR1(:,:,:)) / PD(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_WTH2(:,:,IKB-1)=M3_WTH_WTH2(:,:,IKB)
-M3_WTH_WTH2(:,:,IKE+1)=M3_WTH_WTH2(:,:,IKE)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+M3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD        &
+                   * (1.+0.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)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKB)
+M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKE+1)=M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WTH2',1,ZHOOK_HANDLE)
 END FUNCTION M3_WTH_WTH2
@@ -1011,21 +1051,25 @@ FUNCTION D_M3_WTH_WTH2_O_DDTDZ(D,CSTURB,PM3_WTH_WTH2,PREDTH1,PREDR1,PD,PBLL_O_E,
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_WTH_WTH2_O_DDTDZ
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WTH2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_WTH_WTH2_O_DDTDZ(:,:,:) = (  0.5*CSTURB%XCSHF*PBLL_O_E(:,:,:)*PETHETA(:,:,:)*0.5/CSTURB%XCTD/PD(:,:,:) &
-                                - PM3_WTH_WTH2(:,:,:)/PD(:,:,:)*(1.5+PREDTH1(:,:,:)+PREDR1(:,:,:))  )&
-                             * PBLL_O_E(:,:,:) * PETHETA(:,:,:) * CSTURB%XCTV
-!$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)
+D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = (  0.5*CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PD(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                - PM3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT)/PD(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))  )&
+                             * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
-D_M3_WTH_WTH2_O_DDTDZ(:,:,IKB-1)=D_M3_WTH_WTH2_O_DDTDZ(:,:,IKB)
-D_M3_WTH_WTH2_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_WTH2_O_DDTDZ(:,:,IKE)
+D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WTH2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_WTH_WTH2_O_DDTDZ
@@ -1040,20 +1084,24 @@ FUNCTION M3_WTH_W2TH(D,CSTURB,PREDTH1,PREDR1,PD,PKEFF,PTKE)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_WTH_W2TH
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_W2TH',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-ZWORK1 = MZM(PTKE, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_W2TH(:,:,:) = CSTURB%XCSHF*PKEFF(:,:,:)*1.5/ZWORK1(:,:,:)              &
-  * (1. - 0.5*PREDR1(:,:,:)*(1.+PREDR1(:,:,:))/PD(:,:,:) ) / (1.+PREDTH1(:,:,:))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_WTH_W2TH(:,:,IKB-1)=M3_WTH_W2TH(:,:,IKB)
-M3_WTH_W2TH(:,:,IKE+1)=M3_WTH_W2TH(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+CALL MZM_PHY(D,PTKE,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+M3_WTH_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*1.5/ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)              &
+  * (1. - 0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKB)
+M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKE+1)=M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_W2TH',1,ZHOOK_HANDLE)
 END FUNCTION M3_WTH_W2TH
@@ -1070,24 +1118,28 @@ FUNCTION D_M3_WTH_W2TH_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PBLL_O_E,PETHETA,PKEFF
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_WTH_W2TH_O_DDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK, IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2TH_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-ZWORK1 = MZM(PTKE, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_WTH_W2TH_O_DDTDZ(:,:,:) = &
- - CSTURB%XCSHF*PKEFF(:,:,:)*1.5/ZWORK1(:,:,:)/(1.+PREDTH1(:,:,:))**2 &
- * CSTURB%XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:)  &
- * (1. - 0.5*PREDR1(:,:,:)*(1.+PREDR1(:,:,:))/PD(:,:,:)* &
-   ( 1.+(1.+PREDTH1(:,:,:))*(1.5+PREDR1(:,:,:)+PREDTH1(:,:,:))/PD(:,:,:)) )
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_WTH_W2TH_O_DDTDZ(:,:,IKB-1)=D_M3_WTH_W2TH_O_DDTDZ(:,:,IKB)
-D_M3_WTH_W2TH_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_W2TH_O_DDTDZ(:,:,IKE)
+CALL MZM_PHY(D,PTKE,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = &
+ - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*1.5/ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/(1.+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)  &
+ * (1. - 0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)* &
+   ( 1.+(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.5+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDTH1(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)
+!
+D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2TH_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_WTH_W2TH_O_DDTDZ
@@ -1103,21 +1155,25 @@ FUNCTION M3_WTH_W2R(D,CSTURB,PD,PKEFF,PTKE,PBLL_O_E,PEMOIST,PDTDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_WTH_W2R
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_W2R',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-ZWORK1 = MZM(PTKE, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_W2R(:,:,:) = - CSTURB%XCSHF*PKEFF(:,:,:)*0.75*CSTURB%XCTV*PBLL_O_E(:,:,:) &
-                    /ZWORK1(:,:,:)*PEMOIST(:,:,:)*PDTDZ(:,:,:)/PD(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+CALL MZM_PHY(D,PTKE,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+M3_WTH_W2R(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                    /ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(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)
 !
-M3_WTH_W2R(:,:,IKB-1)=M3_WTH_W2R(:,:,IKB)
-M3_WTH_W2R(:,:,IKE+1)=M3_WTH_W2R(:,:,IKE)
+M3_WTH_W2R(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_W2R(IIB:IIE,IJB:IJE,IKB)
+M3_WTH_W2R(IIB:IIE,IJB:IJE,IKE+1)=M3_WTH_W2R(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_W2R',1,ZHOOK_HANDLE)
 END FUNCTION M3_WTH_W2R
@@ -1134,22 +1190,26 @@ FUNCTION D_M3_WTH_W2R_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PKEFF,PTKE,PBLL_O_E,PEM
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_WTH_W2R_O_DDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2R_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-ZWORK1 = MZM(PTKE, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_WTH_W2R_O_DDTDZ(:,:,:) = - CSTURB%XCSHF*PKEFF(:,:,:)*0.75*CSTURB%XCTV*PBLL_O_E(:,:,:) &
-                               /ZWORK1(:,:,:)*PEMOIST(:,:,:)/PD(:,:,:) &
-                                     * (1. -  PREDTH1(:,:,:)*(1.5+PREDTH1(:,:,:)+PREDR1(:,:,:))/PD(:,:,:))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+CALL MZM_PHY(D,PTKE,ZWORK1)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+                               /ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                     * (1. -  PREDTH1(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))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
-D_M3_WTH_W2R_O_DDTDZ(:,:,IKB-1)=D_M3_WTH_W2R_O_DDTDZ(:,:,IKB)
-D_M3_WTH_W2R_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_W2R_O_DDTDZ(:,:,IKE)
+D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2R_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_WTH_W2R_O_DDTDZ
@@ -1167,22 +1227,28 @@ FUNCTION M3_WTH_WR2(D,CSTURB,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_WTH_WR2
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WR2',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZM(PBETA*PLEPS/(PSQRT_TKE*PTKE), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_WR2(:,:,:) = - CSTURB%XCSHF*PKEFF(:,:,:)*0.25*PBLL_O_E(:,:,:)*CSTURB%XCTV*PEMOIST(:,:,:)**2       &
-                           *ZWORK1(:,:,:)/CSTURB%XCTD*PDTDZ(:,:,:)/PD(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_WTH_WR2(:,:,IKB-1)=M3_WTH_WR2(:,:,IKB)
-M3_WTH_WR2(:,:,IKE+1)=M3_WTH_WR2(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(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)
+M3_WTH_WR2(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2       &
+                           *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*PDTDZ(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)
+!
+M3_WTH_WR2(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_WR2(IIB:IIE,IJB:IJE,IKB)
+M3_WTH_WR2(IIB:IIE,IJB:IJE,IKE+1)=M3_WTH_WR2(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WR2',1,ZHOOK_HANDLE)
 END FUNCTION M3_WTH_WR2
@@ -1201,23 +1267,29 @@ FUNCTION D_M3_WTH_WR2_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PKEFF,PTKE,PSQRT_TKE,PB
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_WTH_WR2_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WR2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZM(PBETA*PLEPS/(PSQRT_TKE*PTKE), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_WTH_WR2_O_DDTDZ(:,:,:) = - CSTURB%XCSHF*PKEFF(:,:,:)*0.25*PBLL_O_E(:,:,:)*CSTURB%XCTV*PEMOIST(:,:,:)**2 &
-                           *ZWORK1(:,:,:)/CSTURB%XCTD/PD(:,:,:)     &
-                           * (1. -  PREDTH1(:,:,:)*(1.5+PREDTH1(:,:,:)+PREDR1(:,:,:))/PD(:,:,:))
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_WTH_WR2_O_DDTDZ(:,:,IKB-1)=D_M3_WTH_WR2_O_DDTDZ(:,:,IKB)
-D_M3_WTH_WR2_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_WR2_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(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)
+D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2 &
+                           *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD/PD(IIB:IIE,IJB:IJE,1:D%NKT)     &
+                           * (1. -  PREDTH1(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))
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WR2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_WTH_WR2_O_DDTDZ
@@ -1234,22 +1306,29 @@ FUNCTION M3_WTH_WTHR(D,CSTURB,PREDR1,PD,PKEFF,PTKE,PSQRT_TKE,PBETA,PLEPS,PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_WTH_WTHR
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WTHR',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-ZWORK1 = MZM(PBETA/PTKE*PSQRT_TKE, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_WTH_WTHR(:,:,:) = CSTURB%XCSHF*PKEFF(:,:,:)*PEMOIST(:,:,:)*ZWORK1(:,:,:) &
-                         *0.5*PLEPS(:,:,:)/CSTURB%XCTD*(1+PREDR1(:,:,:))/PD(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_WTH_WTHR(:,:,IKB-1)=M3_WTH_WTHR(:,:,IKB)
-M3_WTH_WTHR(:,:,IKE+1)=M3_WTH_WTHR(:,:,IKE)
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(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)
+M3_WTH_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
+                         *0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*(1+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)
+!
+M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKB)
+M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKE+1)=M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WTH_WTHR',1,ZHOOK_HANDLE)
 END FUNCTION M3_WTH_WTHR
@@ -1264,20 +1343,24 @@ FUNCTION D_M3_WTH_WTHR_O_DDTDZ(D,CSTURB,PM3_WTH_WTHR,PREDTH1,PREDR1,PD,PBLL_O_E,
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_WTH_WTHR_O_DDTDZ
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WTHR_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_WTH_WTHR_O_DDTDZ(:,:,:) = - PM3_WTH_WTHR(:,:,:) * (1.5+PREDTH1(:,:,:)+PREDR1(:,:,:))&
-                               /PD(:,:,:)*CSTURB%XCTV*PBLL_O_E(:,:,:)*PETHETA(:,:,:)
-!$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)
+D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - PM3_WTH_WTHR(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)*CSTURB%XCTV*PBLL_O_E(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)
 !
-D_M3_WTH_WTHR_O_DDTDZ(:,:,IKB-1)=D_M3_WTH_WTHR_O_DDTDZ(:,:,IKB)
-D_M3_WTH_WTHR_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_WTHR_O_DDTDZ(:,:,IKE)
+D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_WTHR_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_WTH_WTHR_O_DDTDZ
@@ -1293,22 +1376,28 @@ FUNCTION M3_TH2_W2TH(D,CSTURB,PREDTH1,PREDR1,PD,PDTDZ,PLM,PLEPS,PTKE)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_TH2_W2TH
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_W2TH',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((1.-0.5*PREDR1*(1.+PREDR1)/PD)/(1.+PREDTH1)*PDTDZ, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_TH2_W2TH(:,:,:) = - ZWORK1(:,:,:) &
-                       * 1.5*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:)*CSTURB%XCTV
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_TH2_W2TH(:,:,IKB-1)=M3_TH2_W2TH(:,:,IKB)
-M3_TH2_W2TH(:,:,IKE+1)=M3_TH2_W2TH(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))/(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*PDTDZ(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)
+M3_TH2_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = - ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
+                       * 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_TH2_W2TH(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_W2TH(IIB:IIE,IJB:IJE,IKB)
+M3_TH2_W2TH(IIB:IIE,IJB:IJE,IKE+1)=M3_TH2_W2TH(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_W2TH',1,ZHOOK_HANDLE)
 END FUNCTION M3_TH2_W2TH
@@ -1324,36 +1413,44 @@ FUNCTION D_M3_TH2_W2TH_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLM,PLEPS,PTKE,OUSERV)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   LOGICAL,                INTENT(IN) :: OUSERV
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_TH2_W2TH_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_W2TH_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
 
 IF (OUSERV) THEN
-!  D_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM*PLEPS/PTKE*CSTURB%XCTV * MZF(                    &
+!  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM*PLEPS/PTKE*CSTURB%XCTV * MZF(                    &
 !          (1.-0.5*PREDR1*(1.+PREDR1)/PD)*(1.-(1.5+PREDTH1+PREDR1)*(1.+PREDTH1)/PD )  &
 !        / (1.+PREDTH1)**2, D%NKA, D%NKU, D%NKL)
-
-  ZWORK1 = MZF((1.-0.5*PREDR1*(1.+PREDR1)/PD)*(1.-(1.5+PREDTH1+PREDR1)*   &
-             PREDTH1*(1.+PREDTH1)/PD ) / (1.+PREDTH1)**2, D%NKA, D%NKU, D%NKL)
-
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  D_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:)*CSTURB%XCTV * &
-                                 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)
+  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*   &
+             PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))**2
+  !$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)
+  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV * &
+                                 ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
-  ZWORK1 = MZF(1./(1.+PREDTH1)**2, D%NKA, D%NKU, D%NKL)
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-  D_M3_TH2_W2TH_O_DDTDZ(:,:,:) = - 1.5*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:)*CSTURB%XCTV &
-                                 * 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)
+  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = 1./(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))**2
+  !$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)
+  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV &
+                                 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 END IF
 !
-D_M3_TH2_W2TH_O_DDTDZ(:,:,IKB-1)=D_M3_TH2_W2TH_O_DDTDZ(:,:,IKB)
-D_M3_TH2_W2TH_O_DDTDZ(:,:,IKE+1)=D_M3_TH2_W2TH_O_DDTDZ(:,:,IKE)
+D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_W2TH_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_TH2_W2TH_O_DDTDZ
@@ -1367,22 +1464,28 @@ FUNCTION M3_TH2_WTH2(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PSQRT_TKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_TH2_WTH2
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WTH2',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_TH2_WTH2(:,:,:) = PLEPS(:,:,:)*0.5/CSTURB%XCTD/PSQRT_TKE(:,:,:)          &
-                     * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_TH2_WTH2(:,:,IKB-1)=M3_TH2_WTH2(:,:,IKB)
-M3_TH2_WTH2(:,:,IKE+1)=M3_TH2_WTH2(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/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)
+M3_TH2_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PSQRT_TKE(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)
+!
+M3_TH2_WTH2(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_WTH2(IIB:IIE,IJB:IJE,IKB)
+M3_TH2_WTH2(IIB:IIE,IJB:IJE,IKE+1)=M3_TH2_WTH2(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WTH2',1,ZHOOK_HANDLE)
 END FUNCTION M3_TH2_WTH2
@@ -1398,25 +1501,29 @@ FUNCTION D_M3_TH2_WTH2_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_TH2_WTH2_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WTH2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PBLL_O_E*PETHETA* (0.5/PD                                                   &
-             - (1.5+PREDTH1+PREDR1)*(1.+0.5*PREDTH1+1.5*PREDR1+0.5*PREDR1**2)/PD**2 &
-                           ), D%NKA, D%NKU, D%NKL)
- 
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_TH2_WTH2_O_DDTDZ(:,:,:) = PLEPS(:,:,:)*0.5/CSTURB%XCTD/PSQRT_TKE(:,:,:)*CSTURB%XCTV &
-                               * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_TH2_WTH2_O_DDTDZ(:,:,IKB-1)=D_M3_TH2_WTH2_O_DDTDZ(:,:,IKB)
-D_M3_TH2_WTH2_O_DDTDZ(:,:,IKE+1)=D_M3_TH2_WTH2_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)* (0.5/PD(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))*(1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2)
+ !$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)
+D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WTH2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_TH2_WTH2_O_DDTDZ
@@ -1432,22 +1539,28 @@ FUNCTION M3_TH2_W2R(D,CSTURB,PD,PLM,PLEPS,PTKE,PBLL_O_E,PEMOIST,PDTDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_TH2_W2R
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_W2R',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PBLL_O_E*PEMOIST/PD*PDTDZ**2, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_TH2_W2R(:,:,:) = 0.75*CSTURB%XCTV**2*ZWORK1(:,:,:) &
-                    *PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_TH2_W2R(:,:,IKB-1)=M3_TH2_W2R(:,:,IKB)
-M3_TH2_W2R(:,:,IKE+1)=M3_TH2_W2R(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)**2
+!$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)
+M3_TH2_W2R(IIB:IIE,IJB:IJE,1:D%NKT) = 0.75*CSTURB%XCTV**2*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
+                    *PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_TH2_W2R(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_W2R(IIB:IIE,IJB:IJE,IKB)
+M3_TH2_W2R(IIB:IIE,IJB:IJE,IKE+1)=M3_TH2_W2R(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_W2R',1,ZHOOK_HANDLE)
 END FUNCTION M3_TH2_W2R
@@ -1465,22 +1578,28 @@ FUNCTION D_M3_TH2_W2R_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLM,PLEPS,PTKE,PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_TH2_W2R_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_W2R_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 =  MZF(PBLL_O_E*PEMOIST/PD*PDTDZ*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_TH2_W2R_O_DDTDZ(:,:,:) = 0.75*CSTURB%XCTV**2*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:) &
-                              * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_TH2_W2R_O_DDTDZ(:,:,IKB-1)=D_M3_TH2_W2R_O_DDTDZ(:,:,IKB)
-D_M3_TH2_W2R_O_DDTDZ(:,:,IKE+1)=D_M3_TH2_W2R_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(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))
+!$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)
+D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.75*CSTURB%XCTV**2*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(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)
+!
+D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_W2R_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_TH2_W2R_O_DDTDZ
@@ -1495,22 +1614,28 @@ FUNCTION M3_TH2_WR2(D,CSTURB,PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_TH2_WR2
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WR2',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((PBLL_O_E*PEMOIST*PDTDZ)**2/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_TH2_WR2(:,:,:) = 0.25*CSTURB%XCTV**2*ZWORK1(:,:,:)&
-                    *PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_TH2_WR2(:,:,IKB-1)=M3_TH2_WR2(:,:,IKB)
-M3_TH2_WR2(:,:,IKE+1)=M3_TH2_WR2(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT))**2/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)
+M3_TH2_WR2(IIB:IIE,IJB:IJE,1:D%NKT) = 0.25*CSTURB%XCTV**2*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)&
+                    *PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_TH2_WR2(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_WR2(IIB:IIE,IJB:IJE,IKB)
+M3_TH2_WR2(IIB:IIE,IJB:IJE,IKE+1)=M3_TH2_WR2(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WR2',1,ZHOOK_HANDLE)
 END FUNCTION M3_TH2_WR2
@@ -1527,22 +1652,28 @@ FUNCTION D_M3_TH2_WR2_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O_
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_TH2_WR2_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WR2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((PBLL_O_E*PEMOIST)**2*PDTDZ/PD*(2.-PREDTH1*(1.5+PREDTH1+PREDR1)/PD), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_TH2_WR2_O_DDTDZ(:,:,:) = 0.25*CSTURB%XCTV**2*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD &
-                              *  ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_TH2_WR2_O_DDTDZ(:,:,IKB-1)=D_M3_TH2_WR2_O_DDTDZ(:,:,IKB)
-D_M3_TH2_WR2_O_DDTDZ(:,:,IKE+1)=D_M3_TH2_WR2_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT))**2*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(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))
+!$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)
+D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.25*CSTURB%XCTV**2*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
+                              *  ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WR2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_TH2_WR2_O_DDTDZ
@@ -1558,22 +1689,28 @@ FUNCTION M3_TH2_WTHR(D,CSTURB,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_TH2_WTHR
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WTHR',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PBLL_O_E*PEMOIST*PDTDZ*(1.+PREDR1)/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_TH2_WTHR(:,:,:) = - 0.5*CSTURB%XCTV*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD &
-                       *ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_TH2_WTHR(:,:,IKB-1)=M3_TH2_WTHR(:,:,IKB)
-M3_TH2_WTHR(:,:,IKE+1)=M3_TH2_WTHR(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+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)
+M3_TH2_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
+                       *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKB)
+M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKE+1)=M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_TH2_WTHR',1,ZHOOK_HANDLE)
 END FUNCTION M3_TH2_WTHR
@@ -1590,22 +1727,28 @@ FUNCTION D_M3_TH2_WTHR_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PEMOIST
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDTDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_TH2_WTHR_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WTHR_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PBLL_O_E*PEMOIST*(1.+PREDR1)/PD * (1. -PREDTH1*(1.5+PREDTH1+PREDR1)/PD), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_TH2_WTHR_O_DDTDZ(:,:,:) = - 0.5*CSTURB%XCTV*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD &
-                               * ZWORK1(:,:,:) 
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_TH2_WTHR_O_DDTDZ(:,:,IKB-1)=D_M3_TH2_WTHR_O_DDTDZ(:,:,IKB)
-D_M3_TH2_WTHR_O_DDTDZ(:,:,IKE+1)=D_M3_TH2_WTHR_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) * (1. -PREDTH1(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))
+!$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)
+D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_TH2_WTHR_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_TH2_WTHR_O_DDTDZ
@@ -1619,22 +1762,28 @@ FUNCTION M3_THR_WTHR(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PSQRT_TKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_THR_WTHR
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_WTHR',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 =  MZF((1.+PREDTH1)*(1.+PREDR1)/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_THR_WTHR(:,:,:) = 0.5*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD &
-                     * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_THR_WTHR(:,:,IKB-1)=M3_THR_WTHR(:,:,IKB)
-M3_THR_WTHR(:,:,IKE+1)=M3_THR_WTHR(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+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)
+M3_THR_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
+                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_THR_WTHR(IIB:IIE,IJB:IJE,IKB-1)=M3_THR_WTHR(IIB:IIE,IJB:IJE,IKB)
+M3_THR_WTHR(IIB:IIE,IJB:IJE,IKE+1)=M3_THR_WTHR(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_WTHR',1,ZHOOK_HANDLE)
 END FUNCTION M3_THR_WTHR
@@ -1650,22 +1799,28 @@ FUNCTION D_M3_THR_WTHR_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_THR_WTHR_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTHR_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PETHETA*PBLL_O_E/PD*(1.+PREDR1)*(1.-(1.+PREDTH1)*(1.5+PREDTH1+PREDR1)/PD), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_THR_WTHR_O_DDTDZ(:,:,:) = 0.5*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD * CSTURB%XCTV &
-                               * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_THR_WTHR_O_DDTDZ(:,:,IKB-1)=D_M3_THR_WTHR_O_DDTDZ(:,:,IKB)
-D_M3_THR_WTHR_O_DDTDZ(:,:,IKE+1)=D_M3_THR_WTHR_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.+PREDTH1(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))
+!$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)
+D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD * CSTURB%XCTV &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTHR_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_THR_WTHR_O_DDTDZ
@@ -1681,22 +1836,28 @@ FUNCTION M3_THR_WTH2(D,CSTURB,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDRDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDRDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_THR_WTH2
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_WTH2',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((1.+PREDR1)*PBLL_O_E*PETHETA*PDRDZ/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_THR_WTH2(:,:,:) = - 0.25*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD*CSTURB%XCTV &
-                     * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_THR_WTH2(:,:,IKB-1)=M3_THR_WTH2(:,:,IKB)
-M3_THR_WTH2(:,:,IKE+1)=M3_THR_WTH2(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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))*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PDRDZ(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)
+M3_THR_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV &
+                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_THR_WTH2(IIB:IIE,IJB:IJE,IKB-1)=M3_THR_WTH2(IIB:IIE,IJB:IJE,IKB)
+M3_THR_WTH2(IIB:IIE,IJB:IJE,IKE+1)=M3_THR_WTH2(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_WTH2',1,ZHOOK_HANDLE)
 END FUNCTION M3_THR_WTH2
@@ -1713,22 +1874,28 @@ FUNCTION D_M3_THR_WTH2_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDRDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_THR_WTH2_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTH2_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(-(1.+PREDR1)*(PBLL_O_E*PETHETA/PD)**2*PDRDZ*(1.5+PREDTH1+PREDR1), D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_THR_WTH2_O_DDTDZ(:,:,:) = - 0.25*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD*CSTURB%XCTV**2 &
-                               * ZWORK1(:,:,:) 
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_THR_WTH2_O_DDTDZ(:,:,IKB-1)=D_M3_THR_WTH2_O_DDTDZ(:,:,IKB)
-D_M3_THR_WTH2_O_DDTDZ(:,:,IKE+1)=D_M3_THR_WTH2_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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))*(PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT))**2*PDRDZ(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))
+!$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)
+D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV**2 &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTH2_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_THR_WTH2_O_DDTDZ
@@ -1744,24 +1911,29 @@ FUNCTION D_M3_THR_WTH2_O_DDRDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLEPS,PSQRT_TKE,PBLL_O
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PBLL_O_E
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_THR_WTH2_O_DDRDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTH2_O_DDRDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF(PBLL_O_E*PETHETA/PD                                              &
-       *(-(1.+PREDR1)*PREDR1/PD*(1.5+PREDTH1+PREDR1)+(1.+2.*PREDR1)),     &
-       D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_THR_WTH2_O_DDRDZ(:,:,:) = - 0.25*PLEPS(:,:,:)/PSQRT_TKE(:,:,:)/CSTURB%XCTD*CSTURB%XCTV          &
-                               * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_THR_WTH2_O_DDRDZ(:,:,IKB-1)=D_M3_THR_WTH2_O_DDRDZ(:,:,IKB)
-D_M3_THR_WTH2_O_DDRDZ(:,:,IKE+1)=D_M3_THR_WTH2_O_DDRDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)                                              &
+       *(-(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)/PD(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))+(1.+2.*PREDR1(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)
+D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV          &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_WTH2_O_DDRDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_THR_WTH2_O_DDRDZ
@@ -1776,22 +1948,28 @@ FUNCTION M3_THR_W2TH(D,CSTURB,PREDR1,PD,PLM,PLEPS,PTKE,PDRDZ)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDRDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: M3_THR_W2TH
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_W2TH',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 = MZF((1.+PREDR1)*PDRDZ/PD, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-M3_THR_W2TH(:,:,:) = - 0.75*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:) * CSTURB%XCTV      &
-                     * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-M3_THR_W2TH(:,:,IKB-1)=M3_THR_W2TH(:,:,IKB)
-M3_THR_W2TH(:,:,IKE+1)=M3_THR_W2TH(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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))*PDRDZ(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)
+M3_THR_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV      &
+                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+M3_THR_W2TH(IIB:IIE,IJB:IJE,IKB-1)=M3_THR_W2TH(IIB:IIE,IJB:IJE,IKB)
+M3_THR_W2TH(IIB:IIE,IJB:IJE,IKE+1)=M3_THR_W2TH(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_THR_W2TH',1,ZHOOK_HANDLE)
 END FUNCTION M3_THR_W2TH
@@ -1809,22 +1987,28 @@ FUNCTION D_M3_THR_W2TH_O_DDTDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLM,PLEPS,PTKE,PBLL_O_
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDRDZ
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PETHETA
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_THR_W2TH_O_DDTDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_W2TH_O_DDTDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 =  MZF(-PETHETA*PBLL_O_E*(1.+PREDR1)*PDRDZ*(1.5+PREDTH1+PREDR1)/PD**2, D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_THR_W2TH_O_DDTDZ(:,:,:) = - 0.75*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:) * CSTURB%XCTV**2    &
-                               * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_THR_W2TH_O_DDTDZ(:,:,IKB-1)=D_M3_THR_W2TH_O_DDTDZ(:,:,IKB)
-D_M3_THR_W2TH_O_DDTDZ(:,:,IKE+1)=D_M3_THR_W2TH_O_DDTDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+IJB=D%NJBC
+!$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PDRDZ(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
+!$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)
+D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV**2    &
+                               * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_W2TH_O_DDTDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_THR_W2TH_O_DDTDZ
@@ -1839,24 +2023,29 @@ FUNCTION D_M3_THR_W2TH_O_DDRDZ(D,CSTURB,PREDTH1,PREDR1,PD,PLM,PLEPS,PTKE)
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLEPS
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTKE
   REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: D_M3_THR_W2TH_O_DDRDZ
-  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1 ! working array
-  INTEGER :: IKB, IKE, JI,JJ,JK
+  REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: ZWORK1,ZWORK2 ! working array
+  INTEGER :: IKB, IKE, JI,JJ,JK,IIB,IIE,IJB,IJE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_W2TH_O_DDRDZ',0,ZHOOK_HANDLE)
-IKB = 1+JPVEXT_TURB
-IKE = SIZE(PD,3)-JPVEXT_TURB
-
-ZWORK1 =  MZF(-(1.+PREDR1)*PREDR1*(1.5+PREDTH1+PREDR1)/PD**2          &
-        +(1.+2.*PREDR1)/PD,                                    &
-       D%NKA, D%NKU, D%NKL)
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-D_M3_THR_W2TH_O_DDRDZ(:,:,:) = - 0.75*PLM(:,:,:)*PLEPS(:,:,:)/PTKE(:,:,:) * CSTURB%XCTV     &
-                               * ZWORK1(:,:,:)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-!
-D_M3_THR_W2TH_O_DDRDZ(:,:,IKB-1)=D_M3_THR_W2TH_O_DDRDZ(:,:,IKB)
-D_M3_THR_W2TH_O_DDRDZ(:,:,IKE+1)=D_M3_THR_W2TH_O_DDRDZ(:,:,IKE)
+IKB=D%NKTB
+IKE=D%NKTE
+IIE=D%NIEC
+IIB=D%NIBC
+IJE=D%NJEC
+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(:,:,:)
+!$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)
+D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV     &
+                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+!$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
+!
+D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKB)
+D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKE+1)=D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKE)
 !
 IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_THR_W2TH_O_DDRDZ',1,ZHOOK_HANDLE)
 END FUNCTION D_M3_THR_W2TH_O_DDRDZ
-- 
GitLab