From 9ea2495874e1336ab7234b0f9daa57bce7dfa6a4 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 22 Jun 2018 15:55:14 +0200
Subject: [PATCH] Philippe 22/06/2018: retrieval of small optimisations and
 improvements from OPENACC branch

---
 src/MNH/advec_4th_order_aux.f90  |   4 +-
 src/MNH/advec_weno_k_3_aux.f90   |   4 -
 src/MNH/advection_metsv.f90      |   2 +-
 src/MNH/advecuvw_weno_k.f90      |   4 +-
 src/MNH/fast_terms.f90           |  16 ++--
 src/MNH/flat_invz.f90            |  10 ++-
 src/MNH/gradient_m.f90           |   8 +-
 src/MNH/gradient_w.f90           |   6 +-
 src/MNH/phys_paramn.f90          |   2 +-
 src/MNH/ppm.f90                  |  32 ++-----
 src/MNH/prandtl.f90              |  19 ++--
 src/MNH/resolved_cloud.f90       |   1 -
 src/MNH/tke_eps_sources.f90      |  14 ++-
 src/MNH/turb.f90                 | 145 ++++++++++++++++++-------------
 src/MNH/turb_hor_dyn_corr.f90    |  34 ++++----
 src/MNH/turb_ver.f90             |   2 -
 src/MNH/turb_ver_dyn_flux.f90    |  39 ++++-----
 src/MNH/turb_ver_thermo_flux.f90 |   4 +-
 18 files changed, 168 insertions(+), 178 deletions(-)

diff --git a/src/MNH/advec_4th_order_aux.f90 b/src/MNH/advec_4th_order_aux.f90
index 5d98dbed9..6bffab256 100644
--- a/src/MNH/advec_4th_order_aux.f90
+++ b/src/MNH/advec_4th_order_aux.f90
@@ -20,7 +20,7 @@ USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PMEANX, PMEANY ! fluxes
+REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PMEANX, PMEANY ! fluxes
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PFIELDT  ! variable at t
 INTEGER,                INTENT(IN)    :: KGRID    ! C grid localisation
 !
@@ -119,7 +119,7 @@ IMPLICIT NONE
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PMEANX, PMEANY ! fluxes
+REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PMEANX, PMEANY ! fluxes
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PFIELDT  ! variable at t
 INTEGER,                INTENT(IN)    :: KGRID    ! C grid localisation
 !
diff --git a/src/MNH/advec_weno_k_3_aux.f90 b/src/MNH/advec_weno_k_3_aux.f90
index 0115668ae..860a9638f 100644
--- a/src/MNH/advec_weno_k_3_aux.f90
+++ b/src/MNH/advec_weno_k_3_aux.f90
@@ -1493,8 +1493,6 @@ REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)):: ZBNEG1, ZBNEG2, ZBNEG3
 !
 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)):: ZOMP1, ZOMP2, ZOMP3
 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)):: ZOMN1, ZOMN2, ZOMN3
-
-REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZWORK
 !
 ! EPSILON for weno weights calculation
 ! 
@@ -1532,8 +1530,6 @@ ZOMP3  = 0.0
 ZOMN1  = 0.0
 ZOMN2  = 0.0
 ZOMN3  = 0.0 
-
-ZWORK = 0.0
 !
 !-------------------------------------------------------------------------------
 !*       1.1.   Interior Fluxes 
diff --git a/src/MNH/advection_metsv.f90 b/src/MNH/advection_metsv.f90
index 15e583d34..3ec0b2a5c 100644
--- a/src/MNH/advection_metsv.f90
+++ b/src/MNH/advection_metsv.f90
@@ -461,7 +461,7 @@ ZRWCPPM = ZRWCPPM*ZTSTEP_PPM
 !  Clouds    related processes from previous time-step are     taken into account in PRTHS_CLD
 !  Advection related processes from previous time-step will be taken into account in ZRTHS_PPM
 !
-ZRTHS_OTHER = PRTHS - PTHT * PRHODJ / PTSTEP                      
+ZRTHS_OTHER = PRTHS - PTHT * PRHODJ / PTSTEP
 IF (GTKE) ZRTKES_OTHER = PRTKES - PTKET * PRHODJ / PTSTEP                      
 DO JR = 1, KRR
  ZRRS_OTHER(:,:,:,JR) = PRRS(:,:,:,JR) - PRT(:,:,:,JR) * PRHODJ(:,:,:) / PTSTEP
diff --git a/src/MNH/advecuvw_weno_k.f90 b/src/MNH/advecuvw_weno_k.f90
index a28258d06..44c2dd4a6 100644
--- a/src/MNH/advecuvw_weno_k.f90
+++ b/src/MNH/advecuvw_weno_k.f90
@@ -96,7 +96,7 @@ TYPE(HALO2LIST_ll), POINTER :: TZHALO2_UT,TZHALO2_VT,TZHALO2_WT
 TYPE(LIST_ll), POINTER :: TZHALO2_ZMEAN
 INTEGER                     :: IINFO_ll    ! return code of parallel routine
 !
-REAL, DIMENSION(SIZE(PUT,1), SIZE(PUT,2), SIZE(PUT,3)) :: ZMEAN, ZWORK, DYM_ZMEAN
+REAL, DIMENSION(SIZE(PUT,1), SIZE(PUT,2), SIZE(PUT,3)) :: ZMEAN, ZWORK
 !
 INTEGER :: K_SCHEME
 INTEGER :: IKU
@@ -110,8 +110,6 @@ TZHALO2_VT => TPHALO2LIST%NEXT              ! 2nd  add3dfield in model_n
 TZHALO2_WT => TPHALO2LIST%NEXT%NEXT         ! 3rst add3dfield in model_n
 !
 IKU=SIZE(PUT,3)
-ZMEAN = 0.0
-ZWORK=0.0
 !      -------------------------------------------------------
 !
 SELECT CASE(KWENO_ORDER)
diff --git a/src/MNH/fast_terms.f90 b/src/MNH/fast_terms.f90
index 1362c70cd..12bb06a95 100644
--- a/src/MNH/fast_terms.f90
+++ b/src/MNH/fast_terms.f90
@@ -19,11 +19,11 @@ INTERFACE
          !
 INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
 INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
-CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
+CHARACTER(LEN=*),         INTENT(IN)    :: HTURBDIM ! Dimensionality of the
                                                     ! turbulence scheme
-CHARACTER(LEN=4),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
-CHARACTER(LEN=4),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
-CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
+CHARACTER(LEN=*),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
+CHARACTER(LEN=*),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
+CHARACTER(LEN=*),         INTENT(IN)    :: HRAD     ! Radiation scheme name
 LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                     ! Condensation
 REAL,                     INTENT(IN)    :: PTSTEP   ! Time step          
@@ -175,11 +175,11 @@ IMPLICIT NONE
 !
 INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
 INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
-CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
+CHARACTER(LEN=*),         INTENT(IN)    :: HTURBDIM ! Dimensionality of the
                                                     ! turbulence scheme
-CHARACTER(LEN=4),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
-CHARACTER(LEN=4),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
-CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
+CHARACTER(LEN=*),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
+CHARACTER(LEN=*),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
+CHARACTER(LEN=*),         INTENT(IN)    :: HRAD     ! Radiation scheme name
 LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                     ! Condensation
 REAL,                     INTENT(IN)    :: PTSTEP   ! Time step          
diff --git a/src/MNH/flat_invz.f90 b/src/MNH/flat_invz.f90
index 113a7ec4a..e253264a9 100644
--- a/src/MNH/flat_invz.f90
+++ b/src/MNH/flat_invz.f90
@@ -1006,10 +1006,12 @@ CONTAINS
 
   SUBROUTINE FAST_SUBSTITUTION_3D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
        ,PBAND_Y,KIY,KJY,KKU)
-    INTEGER                        :: KIY,KJY,KKU
-    REAL, DIMENSION (KIY*KJY,KKU)  :: PBAND_YR,PBAND_Y,PPBF,PGAM,PAF
-    REAL, DIMENSION (KIY*KJY)      :: PBETX
-    REAL, DIMENSION (KKU)          :: PPCF
+    INTEGER,                      INTENT(IN)  :: KIY,KJY,KKU
+    REAL, DIMENSION (KIY*KJY,KKU),INTENT(OUT) :: PBAND_YR,PGAM
+    REAL, DIMENSION (KIY*KJY,KKU),INTENT(IN)  :: PBAND_Y,PPBF,PAF
+    REAL, DIMENSION (KIY*KJY),    INTENT(OUT) :: PBETX
+    REAL, DIMENSION (KKU),        INTENT(IN)  :: PPCF
+    !
     INTEGER                        :: JK
     !
     !
diff --git a/src/MNH/gradient_m.f90 b/src/MNH/gradient_m.f90
index 7a6abf734..3fed55910 100644
--- a/src/MNH/gradient_m.f90
+++ b/src/MNH/gradient_m.f90
@@ -166,7 +166,7 @@ END MODULE MODI_GRADIENT_M
 !
 !
 USE MODI_SHUMAN
-USE MODD_CONF
+USE MODD_CONF, ONLY:LFLAT
 !
 IMPLICIT NONE
 !
@@ -262,7 +262,7 @@ END FUNCTION GX_M_M
 !*       0.    DECLARATIONS
 !
 !
-USE MODD_CONF
+USE MODD_CONF, ONLY:LFLAT
 USE MODI_SHUMAN
 !
 IMPLICIT NONE
@@ -451,7 +451,7 @@ END FUNCTION GZ_M_M
 !              ------------
 !
 USE MODI_SHUMAN
-USE MODD_CONF
+USE MODD_CONF, ONLY:LFLAT
 USE MODD_PARAMETERS
 !
 IMPLICIT NONE
@@ -604,7 +604,7 @@ END FUNCTION GX_M_U
 !              ------------
 !
 USE MODI_SHUMAN
-USE MODD_CONF
+USE MODD_CONF, ONLY:LFLAT
 USE MODD_PARAMETERS
 !
 IMPLICIT NONE
diff --git a/src/MNH/gradient_w.f90 b/src/MNH/gradient_w.f90
index 80bdd7dc8..a1bfeed20 100644
--- a/src/MNH/gradient_w.f90
+++ b/src/MNH/gradient_w.f90
@@ -40,14 +40,14 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_W_UW ! result UW point
 END FUNCTION GX_W_UW
 !
 !            
-FUNCTION GY_W_VW(KKA,KKU,KL,PA,PDXX,PDZZ,PDZX)      RESULT(PGY_W_VW)
+FUNCTION GY_W_VW(KKA,KKU,KL,PA,PDYY,PDZZ,PDZY)      RESULT(PGY_W_VW)
 !
 INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
 INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the W point
-REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDXX    ! metric coefficient dxx
+REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDYY    ! metric coefficient dyy
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
-REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZX    ! metric coefficient dzx
+REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZY    ! metric coefficient dzy
 !
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_W_VW ! result VW point
 !
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index a50f07b01..bb3e233c2 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -1371,7 +1371,7 @@ END IF
       XTSTEP,TPFILE,                                                        &
       XDXX,XDYY,XDZZ,XDZX,XDZY,XZZ,                                         &
       XDIRCOSXW,XDIRCOSYW,XDIRCOSZW,XCOSSLOPE,XSINSLOPE,                    &
-      XRHODJ,XTHVREF,XRHODREF,                                              &
+      XRHODJ,XTHVREF,                                                       &
       ZSFTH,ZSFRV,ZSFSV,ZSFU,ZSFV,                                          &
       XPABST,XUT,XVT,XWT,XTKET,XSVT,XSRCT,XBL_DEPTH,XSBL_DEPTH,             &
       XCEI,XCEI_MIN,XCEI_MAX,XCOEF_AMPL_SAT,                                &
diff --git a/src/MNH/ppm.f90 b/src/MNH/ppm.f90
index c27623d9a..5ce4d9a1f 100644
--- a/src/MNH/ppm.f90
+++ b/src/MNH/ppm.f90
@@ -232,7 +232,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL0,ZQR0,ZQ60
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
 !BEG JUAN PPM_LL
-INTEGER                          :: ILUOUT,IRESP             ! for prints
 INTEGER                          :: IJS,IJN
 !END JUAN PPM_LL
 !-------------------------------------------------------------------------------
@@ -249,15 +248,6 @@ IJN=IJE
 !*              initialise & update halo & halo2 for PSRC
 !
 CALL GET_HALO(PSRC)
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
 ZFPOS=PSRC
 ZFNEG=PSRC
 !
@@ -347,7 +337,7 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCX(1) == HLBCX(2)
       ZQ6 = 3.0*(ZQL0 - PSRC)
       ZQR = ZQL0 - ZQ6
       ZQL = ZQL0
-   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) 
+   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
       ZQ6 = 3.0*(ZQR0 - PSRC)
       ZQL = ZQR0 - ZQ6
       ZQR = ZQR0
@@ -474,7 +464,7 @@ CASE('OPEN')
 !
    ZQL = ZQL0
    ZQR = ZQR0
-   ZQ6 = ZQ60 
+   ZQ6 = ZQ60
 !
 ! eliminate over and undershoots and create qL and qR as in Lin96
 !
@@ -486,7 +476,7 @@ CASE('OPEN')
       ZQ6 = 3.0*(ZQL0 - PSRC)
       ZQR = ZQL0 - ZQ6
       ZQL = ZQL0
-   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) 
+   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
       ZQ6 = 3.0*(ZQR0 - PSRC)
       ZQL = ZQR0 - ZQ6
       ZQR = ZQR0
@@ -671,7 +661,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL0,ZQR0,ZQ60
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
 !
 !BEG JUAN PPM_LL
-INTEGER                          :: ILUOUT,IRESP             ! for prints
 INTEGER                          :: IIW,IIA
 !END JUAN PPM_LL
 !-------------------------------------------------------------------------------
@@ -688,15 +677,6 @@ CALL GET_HALO(PSRC)
 !-------------------------------------------------------------------------------
 !
 !
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
 ZFPOS=PSRC
 ZFNEG=PSRC
 !
@@ -781,7 +761,7 @@ CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
       ZQ6 = 3.0*(ZQL0 - PSRC)
       ZQR = ZQL0 - ZQ6
       ZQL = ZQL0
-   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) 
+   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
       ZQ6 = 3.0*(ZQR0 - PSRC)
       ZQL = ZQR0 - ZQ6
       ZQR = ZQR0
@@ -903,7 +883,7 @@ CASE('OPEN')
       ZQ6 = 3.0*(ZQL0 - PSRC)
       ZQR = ZQL0 - ZQ6
       ZQL = ZQL0
-   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 ) 
+   ELSEWHERE ( ZQ60*ZDQ > (ZDQ)**2 )
       ZQ6 = 3.0*(ZQR0 - PSRC)
       ZQL = ZQR0 - ZQ6
       ZQR = ZQR0
@@ -1312,7 +1292,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT
 !
 !BEG JUAN PPM_LL
 TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
-INTEGER                          :: ILUOUT,IRESP             ! for prints
 INTEGER                          :: IJS,IJN
 !END JUAN PPM_LL
 !-------------------------------------------------------------------------------
@@ -1577,7 +1556,6 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT
 !BEG JUAN PPM_LL
 TYPE(HALO2LIST_ll), POINTER      :: TZ_PSRC_HALO2_ll         ! halo2 for PSRC
 TYPE(HALO2LIST_ll), POINTER      :: TZ_PHAT_HALO2_ll         ! halo2 for ZPHAT
-INTEGER                          :: ILUOUT,IRESP             ! for prints
 INTEGER                          :: IIW,IIA
 !END JUAN PPM_LL
 !
diff --git a/src/MNH/prandtl.f90 b/src/MNH/prandtl.f90
index 404a45952..13537544f 100644
--- a/src/MNH/prandtl.f90
+++ b/src/MNH/prandtl.f90
@@ -263,7 +263,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PEMOIST ! coefficient E_moist
 !       0.2  declaration of local variables
 !
 REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  &
-                  ZW1, ZW2, ZW3
+                  ZW1, ZW2
 !                                                 working variables
 !                                                     
 INTEGER :: IKB      ! vertical index value for the first inner mass point
@@ -440,19 +440,20 @@ END IF   ! end of the if structure on the turbulence dimensionnality
 !
 !           5. Prandtl numbers for scalars
 !              ---------------------------
-DO JSV=1,ISV
-!
-  IF(HTURBDIM=='1DIM') THEN
+IF(HTURBDIM=='1DIM') THEN
 !        1D case
+  DO JSV=1,ISV
     PRED2THS3(:,:,:,JSV)  = PREDS1(:,:,:,JSV) * PREDTH1(:,:,:)
     IF (KRR /= 0) THEN
       PRED2RS3(:,:,:,JSV)   = PREDR1(:,:,:) *PREDS1(:,:,:,JSV)
     ELSE
       PRED2RS3(:,:,:,JSV)   = 0.
     END IF
+  ENDDO
 !
-  ELSE  IF (L2D) THEN ! 3D case in a 2D model
+ELSE  IF (L2D) THEN ! 3D case in a 2D model
 !
+  DO JSV=1,ISV
     IF (KRR /= 0) THEN
       ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
     ELSE
@@ -473,9 +474,11 @@ DO JSV=1,ISV
     ELSE
       PRED2RS3(:,:,:,JSV) = 0.
     END IF
+  ENDDO
 !
-  ELSE ! 3D case in a 3D model
+ELSE ! 3D case in a 3D model
 !
+  DO JSV=1,ISV
     IF (KRR /= 0) THEN
       ZW1 = MZM(KKA,KKU,KKL, (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
     ELSE
@@ -500,10 +503,10 @@ DO JSV=1,ISV
     ELSE
       PRED2RS3(:,:,:,JSV) = 0.
     END IF
+  ENDDO
 !
-  END IF ! end of HTURBDIM if-block
+END IF ! end of HTURBDIM if-block
 !
-END DO
 !
 !---------------------------------------------------------------------------
 !
diff --git a/src/MNH/resolved_cloud.f90 b/src/MNH/resolved_cloud.f90
index 38e8b8c93..1cd5eb370 100644
--- a/src/MNH/resolved_cloud.f90
+++ b/src/MNH/resolved_cloud.f90
@@ -436,7 +436,6 @@ INTEGER :: IKE           !
 INTEGER :: IKU
 INTEGER :: IINFO_ll      ! return code of parallel routine
 INTEGER :: JK,JI,JL
-INTEGER :: I, J, K
 !
 !
 !
diff --git a/src/MNH/tke_eps_sources.f90 b/src/MNH/tke_eps_sources.f90
index bedd4c18b..883dabdc2 100644
--- a/src/MNH/tke_eps_sources.f90
+++ b/src/MNH/tke_eps_sources.f90
@@ -40,7 +40,7 @@ LOGICAL,                 INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
                                   ! diagnostic fields in the syncronous FM-file
 REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP          ! Dyn. prod. of TKE
 REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTRH
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PTP          ! Ther. prod. of TKE
+REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTP          ! Ther. prod. of TKE
 REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTKES       ! RHOD * Jacobian *
                                                        ! TKE at t+deltat
 REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRTKESM      ! Advection source 
@@ -97,8 +97,6 @@ END MODULE MODI_TKE_EPS_SOURCES
 !!      MXF,MXM.MYF,MYM,MZF,MZM:  Shuman functions (mean operators)
 !!      DZF                    :  Shuman functions (difference operators)     
 !!
-!!      SUBROUTINE TRIDIAG     :  to solve an implicit temporal scheme
-!!      
 !!
 !!    IMPLICIT ARGUMENTS
 !!    ------------------
@@ -192,7 +190,6 @@ USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_SHUMAN 
-USE MODI_TRIDIAG 
 USE MODI_TRIDIAG_TKE
 USE MODI_BUDGET
 USE MODI_LES_MEAN_SUBGRID
@@ -229,8 +226,9 @@ LOGICAL,                 INTENT(IN)   ::  OCLOSE_OUT   ! switch for syncronous
                                                        ! file opening
 LOGICAL,                 INTENT(IN)   ::  OTURB_DIAG   ! switch to write some
                                   ! diagnostic fields in the syncronous FM-file
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP, PTRH          ! Dyn. prod. of TKE
-REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PTP          ! Ther. prod. of TKE
+REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PDP          ! Dyn. prod. of TKE
+REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTRH
+REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PTP          ! Ther. prod. of TKE
 REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTKES       ! RHOD * Jacobian *
                                                        ! TKE at t+deltat
 REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRTHLS       ! Source of Theta_l
@@ -253,8 +251,8 @@ REAL, DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3))::         &
        ZFLX,     & ! horizontal or vertical flux of the treated variable
        ZSOURCE,  & ! source of evolution for the treated variable
        ZKEFF       ! effectif diffusion coeff = LT * SQRT( TKE )
-LOGICAL,DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3)) :: GTKENEG
-                   ! 3D mask .T. if TKE < XTKEMIN
+!LOGICAL,DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3)) :: GTKENEG
+!                   ! 3D mask .T. if TKE < XTKEMIN
 INTEGER             :: IIB,IIE,IJB,IJE,IKB,IKE
                                     ! Index values for the Beginning and End
                                     ! mass points of the domain 
diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90
index 5b4090e65..5e304989b 100644
--- a/src/MNH/turb.f90
+++ b/src/MNH/turb.f90
@@ -15,7 +15,7 @@ INTERFACE
                 HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HCLOUD,PIMPL,      &
                 PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,           &
                 PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
-                PRHODJ,PTHVREF,PRHODREF,                              &
+                PRHODJ,PTHVREF,                                       &
                 PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
                 PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
                 PBL_DEPTH, PSBL_DEPTH,                                &
@@ -70,8 +70,6 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE       ! sinus of the angle
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODJ    ! dry density * Grid size
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PTHVREF   ! Virtual Potential
                                         ! Temperature of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODREF  ! dry density of the 
-                                        ! reference state
 !
 REAL, DIMENSION(:,:),   INTENT(IN)      ::  PSFTH,PSFRV,   &
 ! normal surface fluxes of theta and Rv 
@@ -121,11 +119,11 @@ REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PFLXZTHVMF
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTH       ! heat flux
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWRC       ! cloud water flux
 REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PDYP  ! Dynamical production of TKE
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PTHP  ! Thermal production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PTR   ! Transport production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PDISS ! Dissipation of TKE
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PLEM  ! Mixing length                   
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDYP  ! Dynamical production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHP  ! Thermal production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTR   ! Transport production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDISS ! Dissipation of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PLEM  ! Mixing length
 
 !
 !-------------------------------------------------------------------------------
@@ -142,7 +140,7 @@ END MODULE MODI_TURB
                 HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HCLOUD,PIMPL,      &
                 PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,           &
                 PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE,    &
-                PRHODJ,PTHVREF,PRHODREF,                              &
+                PRHODJ,PTHVREF,                                       &
                 PSFTH,PSFRV,PSFSV,PSFU,PSFV,                          &
                 PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT,                  &
                 PBL_DEPTH,PSBL_DEPTH,                                 &
@@ -432,13 +430,11 @@ REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSINSLOPE       ! sinus of the angle
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODJ    ! dry density * Grid size
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PTHVREF   ! Virtual Potential
                                         ! Temperature of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PRHODREF  ! dry density of the 
-                                        ! reference state
 !
 REAL, DIMENSION(:,:),   INTENT(IN)      ::  PSFTH,PSFRV,   &
 ! normal surface fluxes of theta and Rv 
                                             PSFU,PSFV
-! normal surface fluxes of (u,v) parallel to the orography 
+! normal surface fluxes of (u,v) parallel to the orography
 REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PSFSV
 ! normal surface fluxes of Scalar var. 
 !
@@ -483,11 +479,11 @@ REAL, DIMENSION(:,:,:), INTENT(IN)      ::  PFLXZTHVMF
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWTH       ! heat flux
 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PWRC       ! cloud water flux
 REAL, DIMENSION(:,:,:,:),INTENT(OUT) :: PWSV       ! scalar flux
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PDYP  ! Dynamical production of TKE
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PTHP  ! Thermal production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PTR   ! Transport production of TKE
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PDISS ! Dissipation of TKE
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PLEM  ! Mixing length                   
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDYP  ! Dynamical production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHP  ! Thermal production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTR   ! Transport production of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PDISS ! Dissipation of TKE
+REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PLEM  ! Mixing length
 !
 !
 !-------------------------------------------------------------------------------
@@ -619,7 +615,7 @@ ZRM(:,:,:,:) = PRT(:,:,:,:)
 !
 !*      2.1 Cph at t
 !
-ZCP=XCPD
+ZCP(:,:,:)=XCPD
 !
 IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PRT(:,:,:,1)
 DO JRR = 2,1+KRRL                          ! loop on the liquid components  
@@ -714,7 +710,7 @@ END IF              ! loop end on KRRL >= 1
 ! computes conservative variables
 !
 IF ( KRRL >= 1 ) THEN
-  IF ( KRRI >= 1 ) THEN 
+  IF ( KRRI >= 1 ) THEN
     ! Rnp at t
     PRT(:,:,:,1)  = PRT(:,:,:,1)  + PRT(:,:,:,2)  + PRT(:,:,:,4)
     PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2) + PRRS(:,:,:,4)
@@ -801,7 +797,7 @@ IF (KMODEL_CL==KMI .AND. HTURBLEN_CL/='NONE') CALL CLOUD_MODIF_LM
 !*      3.6 Dissipative length
 !           ------------------
 !
-ZLEPS=ZLM
+ZLEPS(:,:,:)=ZLM(:,:,:)
 !
 !*      3.7 Correction in the Surface Boundary Layer (Redelsperger 2001)
 !           ----------------------------------------
@@ -876,24 +872,35 @@ ZMTH2 = 0.     ! w'th'2
 ZMR2  = 0.     ! w'r'2
 ZMTHR = 0.     ! w'th'r'
 
-IF (HTOM=='TM06') CALL TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
-!
-ZFWTH = -GZ_M_W(KKA,KKU,KKL,ZMWTH,PDZZ)    ! -d(w'2th' )/dz
-ZFWR  = -GZ_M_W(KKA,KKU,KKL,ZMWR, PDZZ)    ! -d(w'2r'  )/dz
-ZFTH2 = -GZ_W_M(KKA,KKU,KKL,ZMTH2,PDZZ)    ! -d(w'th'2 )/dz
-ZFR2  = -GZ_W_M(KKA,KKU,KKL,ZMR2, PDZZ)    ! -d(w'r'2  )/dz
-ZFTHR = -GZ_W_M(KKA,KKU,KKL,ZMTHR,PDZZ)    ! -d(w'th'r')/dz
-!
-ZFWTH(:,:,IKTE:) = 0.
-ZFWTH(:,:,:IKTB) = 0.
-ZFWR (:,:,IKTE:) = 0.
-ZFWR (:,:,:IKTB) = 0.
-ZFTH2(:,:,IKTE:) = 0.
-ZFTH2(:,:,:IKTB) = 0.
-ZFR2 (:,:,IKTE:) = 0.
-ZFR2 (:,:,:IKTB) = 0.
-ZFTHR(:,:,IKTE:) = 0.
-ZFTHR(:,:,:IKTB) = 0.
+IF (HTOM=='TM06') THEN
+  CALL TM06(KKA,KKU,KKL,PTHVREF,PBL_DEPTH,PZZ,PSFTH,ZMWTH,ZMTH2)
+!
+  ZFWTH = -GZ_M_W(KKA,KKU,KKL,ZMWTH,PDZZ)    ! -d(w'2th' )/dz
+  !ZFWR  = -GZ_M_W(KKA,KKU,KKL,ZMWR, PDZZ)    ! -d(w'2r'  )/dz
+  ZFTH2 = -GZ_W_M(KKA,KKU,KKL,ZMTH2,PDZZ)    ! -d(w'th'2 )/dz
+  !ZFR2  = -GZ_W_M(KKA,KKU,KKL,ZMR2, PDZZ)    ! -d(w'r'2  )/dz
+  !ZFTHR = -GZ_W_M(KKA,KKU,KKL,ZMTHR,PDZZ)    ! -d(w'th'r')/dz
+!
+  ZFWTH(:,:,IKTE:) = 0.
+  ZFWTH(:,:,:IKTB) = 0.
+  !ZFWR (:,:,IKTE:) = 0.
+  !ZFWR (:,:,:IKTB) = 0.
+  ZFWR  = 0.
+  ZFTH2(:,:,IKTE:) = 0.
+  ZFTH2(:,:,:IKTB) = 0.
+  !ZFR2 (:,:,IKTE:) = 0.
+  !ZFR2 (:,:,:IKTB) = 0.
+  ZFR2  = 0.
+  !ZFTHR(:,:,IKTE:) = 0.
+  !ZFTHR(:,:,:IKTB) = 0.
+  ZFTHR = 0.
+ELSE
+  ZFWTH = 0.
+  ZFWR  = 0.
+  ZFTH2 = 0.
+  ZFR2  = 0.
+  ZFTHR = 0.
+ENDIF
 !
 !----------------------------------------------------------------------------
 !
@@ -1171,8 +1178,8 @@ IF (LLES_CALL) THEN
 !
   IF (HTURBDIM=="1DIM") THEN
     CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_U2)
-    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_V2)
-    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_W2)
+    X_LES_SUBGRID_V2 = X_LES_SUBGRID_U2
+    X_LES_SUBGRID_W2 = X_LES_SUBGRID_U2
     CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(KKA,KKU,KKL,&
                & GZ_M_W(KKA,KKU,KKL,PTHLT,PDZZ)),X_LES_RES_ddz_Thl_SBG_W2)
     IF (KRR>=1) &
@@ -1308,7 +1315,7 @@ IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments 
 !
-REAL                                  :: PALP,PBETA,PGAM,PLTT,PC
+REAL,                   INTENT(IN)    :: PALP,PBETA,PGAM,PLTT,PC
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PT,PEXN,PCP
 !
 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLOCPEXN
@@ -1463,7 +1470,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLM
 !*       0.2   Declarations of local variables
 !
 REAL                :: ZD           ! distance to the surface
-REAL, DIMENSION(:,:), ALLOCATABLE  ::   ZWORK2D
+REAL                :: ZVAR         ! Intermediary variable
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) ::   ZWORK2D
 !
 REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) ::     &
             ZDTHLDZ,ZDRTDZ,     &!dtheta_l/dz, drt_dz used for computing the stablity
@@ -1488,29 +1496,44 @@ IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
 END IF
 !   compute a mixing length limited by the stability
 !
-ALLOCATE(ZWORK2D(SIZE(PUT,1),SIZE(PUT,2)))
-!
 ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
 ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
 !
-DO JK = IKTB+1,IKTE-1
-      ZDTHLDZ(:,:,JK)= 0.5*((PTHLT(:,:,JK+KKL)-PTHLT(:,:,JK))/PDZZ(:,:,JK+KKL)+      &
-                         (PTHLT(:,:,JK)-PTHLT(:,:,JK-KKL))/PDZZ(:,:,JK))
 ! For dry simulations
-IF (KRR>0) THEN    
-        ZDRTDZ(:,:,JK)= 0.5*((PRT(:,:,JK+KKL,1)-PRT(:,:,JK,1))/PDZZ(:,:,JK+KKL)+      &
-                         (PRT(:,:,JK,1)-PRT(:,:,JK-KKL,1))/PDZZ(:,:,JK))
+IF (KRR>0) THEN
+  DO JK = IKTB+1,IKTE-1
+    DO JJ=1,SIZE(PUT,2)
+      DO JI=1,SIZE(PUT,1)
+        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
+        ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK    ,1))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PRT(JI,JJ,JK    ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK    ))
+        ZVAR=XG/PTHVREF(JI,JJ,JK)*                                                  &
+             (ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK))
+        !
+        IF (ZVAR>0.) THEN
+          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
+                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
+        END IF
+      END DO
+    END DO
+  END DO
 ELSE
-       ZDRTDZ(:,:,JK)=0
-ENDIF
-       ZWORK2D(:,:)=XG/PTHVREF(:,:,JK)*                                           &
-            (ZETHETA(:,:,JK)*ZDTHLDZ(:,:,JK)+ZEMOIST(:,:,JK)*ZDRTDZ(:,:,JK))
-  !
-  WHERE(ZWORK2D(:,:)>0.)
-    PLM(:,:,JK)=MAX(XMNH_EPSILON,MIN(PLM(:,:,JK),                &
-                    0.76* SQRT(PTKET(:,:,JK)/ZWORK2D(:,:))))
-  END WHERE
-END DO
+  DO JK = IKTB+1,IKTE-1
+    DO JJ=1,SIZE(PUT,2)
+      DO JI=1,SIZE(PUT,1)
+        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
+                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
+        ZVAR=XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
+        !
+        IF (ZVAR>0.) THEN
+          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
+                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
+        END IF
+      END DO
+    END DO
+  END DO
+END IF
 !  special case near the surface 
 ZDTHLDZ(:,:,IKB)=(PTHLT(:,:,IKB+KKL)-PTHLT(:,:,IKB))/PDZZ(:,:,IKB+KKL)
 ! For dry simulations
@@ -1527,8 +1550,6 @@ WHERE(ZWORK2D(:,:)>0.)
                     0.76* SQRT(PTKET(:,:,IKB)/ZWORK2D(:,:))))
 END WHERE
 !
-DEALLOCATE(ZWORK2D)
-!
 !  mixing length limited by the distance normal to the surface (with the same factor as for BL89)
 !
 IF (.NOT. ORMC01) THEN
diff --git a/src/MNH/turb_hor_dyn_corr.f90 b/src/MNH/turb_hor_dyn_corr.f90
index 9409dd225..930b8fb69 100644
--- a/src/MNH/turb_hor_dyn_corr.f90
+++ b/src/MNH/turb_hor_dyn_corr.f90
@@ -351,8 +351,8 @@ ZDW_DZ(:,:,:)=-ZDU_DX(:,:,:)-ZDV_DY(:,:,:)
 !
 !* computation 
 !
-ZFLX(:,:,IKB:IKB)   = (2./3.) * PTKEM(:,:,IKB:IKB)                           &
-  - XCMFS * PK(:,:,IKB:IKB) * 2. * ZDU_DX(:,:,:)
+ZFLX(:,:,IKB)   = (2./3.) * PTKEM(:,:,IKB)                           &
+  - XCMFS * PK(:,:,IKB) * 2. * ZDU_DX(:,:,1)
 
 
 !!  &  to be tested later
@@ -360,17 +360,16 @@ ZFLX(:,:,IKB:IKB)   = (2./3.) * PTKEM(:,:,IKB:IKB)                           &
 !!   (-2./3.) * PTP(:,:,IKB:IKB)
 !
 ! extrapolates this flux under the ground with the surface flux
-ZFLX(:,:,IKB-1) =                                                          &
+ZFLX(:,:,IKB-1) =                                                            &
         PTAU11M(:,:) * PCOSSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2                 &
   -2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:)       &
   +     PTAU22M(:,:) * PSINSLOPE(:,:)**2                                     &
   +     PTAU33M(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2                 &
   +2. * PCDUEFF(:,:) *      (                                                &
       PVSLOPEM(:,:) * PCOSSLOPE(:,:)    * PSINSLOPE(:,:) * ZDIRSINZW(:,:)    &
-    - PUSLOPEM(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)    &
-                            ) 
+    - PUSLOPEM(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)    )
 ! 
-ZFLX(:,:,IKB-1:IKB-1) = 2. * ZFLX(:,:,IKB-1:IKB-1) -  ZFLX(:,:,IKB:IKB)
+ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) -  ZFLX(:,:,IKB)
 !
 CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll)
 IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN
@@ -447,25 +446,24 @@ END IF
 !
 ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE) 
 !
-ZFLX(:,:,IKB:IKB)   = (2./3.) * PTKEM(:,:,IKB:IKB)                           &
-  - XCMFS * PK(:,:,IKB:IKB) * 2. * ZDV_DY(:,:,:)
+ZFLX(:,:,IKB)   = (2./3.) * PTKEM(:,:,IKB)                           &
+  - XCMFS * PK(:,:,IKB) * 2. * ZDV_DY(:,:,1)
 
 !!           & to be tested
 !! + XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) *         &
 !!   (-2./3.) * PTP(:,:,IKB:IKB)
 !
 ! extrapolates this flux under the ground with the surface flux
-ZFLX(:,:,IKB-1) =                                                          &
+ZFLX(:,:,IKB-1) =                                                            &
         PTAU11M(:,:) * PSINSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2                 &         
   +2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:)       &
   +     PTAU22M(:,:) * PCOSSLOPE(:,:)**2                                     &
   +     PTAU33M(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2                 &
   -2. * PCDUEFF(:,:)*       (                                                &
       PUSLOPEM(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)    &
-    + PVSLOPEM(:,:) * PCOSSLOPE(:,:)    * PSINSLOPE(:,:) * ZDIRSINZW(:,:)    &
-                            ) 
+    + PVSLOPEM(:,:) * PCOSSLOPE(:,:)    * PSINSLOPE(:,:) * ZDIRSINZW(:,:)    )
 ! 
-ZFLX(:,:,IKB-1:IKB-1) = 2. * ZFLX(:,:,IKB-1:IKB-1) -  ZFLX(:,:,IKB:IKB)
+ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) -  ZFLX(:,:,IKB)
 !
 CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll)
 !
@@ -541,16 +539,16 @@ END IF
 !
 ZFLX(:,:,IKE+1)= ZFLX(:,:,IKE)
 !
-ZFLX(:,:,IKB:IKB)   = (2./3.) * PTKEM(:,:,IKB:IKB)                           &
-  - XCMFS * PK(:,:,IKB:IKB) * 2. * ZDW_DZ(:,:,:)
+ZFLX(:,:,IKB)   = (2./3.) * PTKEM(:,:,IKB)                           &
+  - XCMFS * PK(:,:,IKB) * 2. * ZDW_DZ(:,:,1)
 
-!!             &  to be tested
-!!   - 2.* XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) *             &
-!!  (-2./3.) * PTP(:,:,IKB:IKB)
+!             &  to be tested
+!   - 2.* XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) *             &
+!  (-2./3.) * PTP(:,:,IKB:IKB)
 !
 ! extrapolates this flux under the ground with the surface flux
 ZFLX(:,:,IKB-1) =                                                     &
-        PTAU11M(:,:) * ZDIRSINZW(:,:)**2                                &         
+        PTAU11M(:,:) * ZDIRSINZW(:,:)**2                                &
   +     PTAU33M(:,:) * PDIRCOSZW(:,:)**2                                &
   +2. * PCDUEFF(:,:)* PUSLOPEM(:,:)  * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) 
   ! 
diff --git a/src/MNH/turb_ver.f90 b/src/MNH/turb_ver.f90
index 51b3acd23..ecf077cbb 100644
--- a/src/MNH/turb_ver.f90
+++ b/src/MNH/turb_ver.f90
@@ -517,8 +517,6 @@ ALLOCATE ( &
 !
 !*       1.   PRELIMINARIES
 !             -------------
-PTP (:,:,:) = 0.
-PDP (:,:,:) = 0.
 !
 IKB=KKA+JPVEXT_TURB*KKL
 IKE=KKU-JPVEXT_TURB*KKL
diff --git a/src/MNH/turb_ver_dyn_flux.f90 b/src/MNH/turb_ver_dyn_flux.f90
index 3f39c0030..61ed8f8be 100644
--- a/src/MNH/turb_ver_dyn_flux.f90
+++ b/src/MNH/turb_ver_dyn_flux.f90
@@ -75,8 +75,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PWV          ! momentum flux v'w'
 REAL, DIMENSION(:,:,:), INTENT(INOUT)   ::  PRUS, PRVS, PRWS
                             ! cumulated sources for the prognostic variables
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT)::  PDP,PTP   ! Dynamic and thermal
-                                                   ! TKE production terms
+REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PDP          ! Dynamic TKE production term
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTP          ! Thermal TKE production term
 !
 !
 !
@@ -365,8 +365,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PWV          ! momentum flux v'w'
 REAL, DIMENSION(:,:,:), INTENT(INOUT)   ::  PRUS, PRVS, PRWS
                             ! cumulated sources for the prognostic variables
 !
-REAL, DIMENSION(:,:,:), INTENT(INOUT)::  PDP,PTP   ! Dynamic and thermal
-                                                   ! TKE production terms
+REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PDP          ! Dynamic TKE production term
+REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTP          ! Thermal TKE production term
 !
 !
 !
@@ -419,8 +419,7 @@ IKTE=IKT-JPVEXT_TURB
 
 
 !
-ZSOURCE = 0.
-ZFLXZ   = 0.
+ZSOURCE(:,:,:) = 0.
 !
 ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
 !  compute the coefficients for the uncentred gradient computation near the 
@@ -463,20 +462,20 @@ ZCOEFS(:,:,1)=  ZCOEFFLXU(:,:,1) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:)  &
 ZCOEFS(:,:,1:1)=MXM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
 !
 ! compute the explicit tangential flux at the W point
-ZSOURCE(:,:,IKB)     =                                                    &
-    PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:)         &
-   -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:)                          &
+ZSOURCE(:,:,IKB)     =                                              &
+    PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+   -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:)                  &
    -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)  
 !
 ! add the vertical part or the surface flux at the U,W vorticity point
 
-ZSOURCE(:,:,IKB:IKB) =                                      &
-  (   MXM( ZSOURCE(:,:,IKB:IKB)   / PDZZ(:,:,IKB:IKB) )     &
-   +  MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB)      &
+ZSOURCE(:,:,IKB:IKB) =                                  &
+  (   MXM( ZSOURCE(:,:,IKB:IKB)   / PDZZ(:,:,IKB:IKB) ) &
+   +  MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB)       &
            *ZUSLOPEM(:,:,1:1)                           &
-          -ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB)      &
-           *ZVSLOPEM(:,:,1:1)                      )     &
-   -  ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL            &
+          -ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB)       &
+           *ZVSLOPEM(:,:,1:1)                      )    &
+   -  ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL        &
   ) * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
 !
 ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
@@ -620,7 +619,7 @@ END IF
 ! Preparation of the arguments for TRIDIAG_WIND
 !!
 ZA(:,:,:)    = - PTSTEP * XCMFS *                              &
-              MYM( ZKEFF ) * MYM(MZM(KKA,KKU,KKL, PRHODJ )) / &
+              MYM( ZKEFF ) * MYM(MZM(KKA,KKU,KKL, PRHODJ )) /  &
               MYM( PDZZ )**2
 !
 !
@@ -648,11 +647,11 @@ ZSOURCE(:,:,IKB)       =                                                  &
 !
 ! add the vertical part or the surface flux at the V,W vorticity point
 ZSOURCE(:,:,IKB:IKB) =                                      &
-  (   MYM( ZSOURCE(:,:,IKB:IKB)   / PDZZ(:,:,IKB:IKB) )         &
+  (   MYM( ZSOURCE(:,:,IKB:IKB)   / PDZZ(:,:,IKB:IKB) )     &
    +  MYM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB)           &
-          *ZUSLOPEM(:,:,1:1)                            &
+          *ZUSLOPEM(:,:,1:1)                                &
           +ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB)           &
-          *ZVSLOPEM(:,:,1:1)                      )     &
+          *ZVSLOPEM(:,:,1:1)                      )         &
    - ZCOEFS(:,:,1:1) * PVM(:,:,IKB:IKB) * PIMPL             &
   ) * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
 !
@@ -760,7 +759,7 @@ IF(HTURBDIM=='3DIM') THEN
                 /(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB  ))           &
             )                                                              &
           * PDZY(:,:,IKB+KKL:IKB+KKL)                                      &
-       ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB)))                     &
+       ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB)))             &
                             )
   !
     PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
diff --git a/src/MNH/turb_ver_thermo_flux.f90 b/src/MNH/turb_ver_thermo_flux.f90
index 7e5206e0d..0915d2544 100644
--- a/src/MNH/turb_ver_thermo_flux.f90
+++ b/src/MNH/turb_ver_thermo_flux.f90
@@ -106,7 +106,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS       ! cumulated source for rt
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTHLP      ! guess of thl at t+ deltat
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRP        ! guess of r at t+ deltat
 !
-REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PTP       ! Dynamic and thermal
+REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTP       ! Dynamic and thermal
                                                      ! TKE production terms
 !
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PWTH       ! heat flux
@@ -438,7 +438,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS       ! cumulated source for rt
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTHLP      ! guess of thl at t+ deltat
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRP        ! guess of r at t+ deltat
 !
-REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PTP       ! Dynamic and thermal
+REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PTP       ! Dynamic and thermal
                                                      ! TKE production terms
 !
 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PWTH       ! heat flux
-- 
GitLab