Index: /LMDZ5/trunk/libf/phylmd/coef_diff_turb_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 1761)
@@ -158,10 +158,10 @@
          
 !   iflag_pbl peut etre utilise comme longuer de melange
-       IF (iflag_pbl.GE.18) THEN
+       IF (iflag_pbl.GE.31) THEN
           CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
                yzlev,yzlay,yu,yv,yteta, &
                ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
                iflag_pbl)
-       ELSE
+       ELSE IF (iflag_pbl<20) THEN
           CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
                yzlev,yzlay,yu,yv,yteta, &
Index: /LMDZ5/trunk/libf/phylmd/ener_conserv.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/ener_conserv.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/ener_conserv.F90	(revision 1761)
@@ -1,4 +1,4 @@
 subroutine ener_conserv(klon,klev,pdtphys, &
- &    puo,pvo,pto,pqo,pun,pvn,ptn,pqn,masse,exner,d_t_ec)
+ &    puo,pvo,pto,pqo,pun,pvn,ptn,pqn,dtke,masse,exner,d_t_ec)
 
 !=============================================================
@@ -19,6 +19,6 @@
 
 ! From module
-USE phys_local_var_mod, ONLY : d_u_vdf,d_v_vdf,d_t_vdf,d_u_ajs,d_v_ajs,d_t_ajs,d_u_con,d_v_con,d_t_con
-USE phys_output_var_mod, ONLY : bils_ec,bils_kinetic,bils_enthalp,bils_latent
+USE phys_local_var_mod, ONLY : d_u_vdf,d_v_vdf,d_t_vdf,d_u_ajs,d_v_ajs,d_t_ajs,d_u_con,d_v_con,d_t_con,d_t_diss
+USE phys_output_var_mod, ONLY : bils_ec,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss
 
 IMPLICIT none
@@ -26,4 +26,5 @@
 #include "YOETHF.h"
 #include "clesphys.h"
+#include "compbl.h"
 
 ! Arguments
@@ -33,4 +34,5 @@
 REAL, DIMENSION(klon,klev),INTENT(IN) :: pun,pvn,ptn,pqn
 REAL, DIMENSION(klon,klev),INTENT(IN) :: masse,exner
+REAL, DIMENSION(klon,klev+1),INTENT(IN) :: dtke
 REAL, DIMENSION(klon,klev),INTENT(OUT) :: d_t_ec
       integer k,i
@@ -65,7 +67,14 @@
 
    IF (iflag_ener_conserv<=2) THEN
-      d_t(:,:)=d_t_vdf(:,:)+d_t_ajs(:,:)   ! d_t_ajs = adjust + thermals
-      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)
-      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)
+!     print*,'ener_conserv pbl=',iflag_pbl
+      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN !d_t_diss accounts for conserv
+         d_t(:,:)=d_t_ajs(:,:)   ! d_t_ajs = adjust + thermals
+         d_u(:,:)=d_u_ajs(:,:)+d_u_con(:,:)
+         d_v(:,:)=d_v_ajs(:,:)+d_v_con(:,:)
+      ELSE
+         d_t(:,:)=d_t_vdf(:,:)+d_t_ajs(:,:)   ! d_t_ajs = adjust + thermals
+         d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)
+         d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)
+      ENDIF
    ELSEIF (iflag_ener_conserv==101) THEN
       d_t(:,:)=0.
@@ -89,6 +98,11 @@
       zv(:,:)=pvo(:,:)
    else
-      zu(:,:)=puo(:,:)+0.5*d_u(:,:)
-      zv(:,:)=pvo(:,:)+0.5*d_v(:,:)
+      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN
+         zu(:,:)=puo(:,:)+d_u_vdf(:,:)+0.5*d_u(:,:)
+         zv(:,:)=pvo(:,:)+d_v_vdf(:,:)+0.5*d_v(:,:)
+      ELSE
+         zu(:,:)=puo(:,:)+0.5*d_u(:,:)
+         zv(:,:)=pvo(:,:)+0.5*d_v(:,:)
+      ENDIF
    endif
 
@@ -130,4 +144,6 @@
 
       bils_ec(:)=0.
+      bils_tke(:)=0.
+      bils_diss(:)=0.
       bils_kinetic(:)=0.
       bils_enthalp(:)=0.
@@ -135,4 +151,6 @@
       DO k=1,klev
         bils_ec(:)=bils_ec(:)-d_t_ec(:,k)*masse(:,k)
+        bils_tke(:)=bils_tke(:)+0.5*(dtke(:,k)+dtke(:,k+1))*masse(:,k)
+        bils_diss(:)=bils_diss(:)-d_t_diss(:,k)*masse(:,k)
         bils_kinetic(:)=bils_kinetic(:)+masse(:,k)* &
      &           (pun(:,k)*pun(:,k)+pvn(:,k)*pvn(:,k) &
@@ -144,9 +162,9 @@
       ENDDO
       bils_ec(:)=rcpd*bils_ec(:)/pdtphys
+      bils_tke(:)=bils_tke(:)/pdtphys
+      bils_diss(:)=rcpd*bils_diss(:)/pdtphys
       bils_kinetic(:)= 0.5*bils_kinetic(:)/pdtphys
       bils_enthalp(:)=rcpd*bils_enthalp(:)/pdtphys
       bils_latent(:)=rlvtt*bils_latent(:)/pdtphys
-
-
 RETURN
 
Index: /LMDZ5/trunk/libf/phylmd/indicesol.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/indicesol.h	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/indicesol.h	(revision 1761)
@@ -21,4 +21,7 @@
       PARAMETER (is_lic=2) ! glacier continental
 !
+      INTEGER is_ave
+      PARAMETER (is_ave=nbsrf+1) ! glacier continental
+!
       REAL epsfra
       PARAMETER (epsfra=1.0E-05)
Index: /LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 1761)
@@ -176,5 +176,5 @@
        alb1_m,    alb2_m,    zxsens,   zxevap,        &
        zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
-       d_t,       d_q,       d_u,      d_v,           & 
+       d_t,       d_q,       d_u,      d_v, d_t_diss, & 
        zcoefh,    zcoefm,    slab_wfbils,             &
        qsol_d,    zq2m,      s_pblh,   s_plcl,        &
@@ -259,4 +259,6 @@
     INCLUDE "YOETHF.h"
     INCLUDE "temps.h"
+!****************************************************************************************
+! Declarations specifiques pour le 1D. A reprendre 
 ! Input variables
 !****************************************************************************************
@@ -291,5 +293,5 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
-    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke 
+    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke 
 
 ! Output variables
@@ -310,9 +312,10 @@
     REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature 
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t_diss       ! change in temperature 
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_q        ! change in water vapour
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
-    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh     ! coef for turbulent diffusion of T and Q, mean for each grid point
-    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefm     ! coef for turbulent diffusion of U and V (?), mean for each grid point
+    REAL, DIMENSION(klon, klev,nbsrf+1),  INTENT(OUT)       :: zcoefh     ! coef for turbulent diffusion of T and Q, mean for each grid point
+    REAL, DIMENSION(klon, klev,nbsrf+1),  INTENT(OUT)       :: zcoefm     ! coef for turbulent diffusion of U and V (?), mean for each grid point
 
 ! Output only for diagnostics
@@ -424,9 +427,9 @@
     REAL, DIMENSION(klon)              :: ztsol
     REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
-    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q
+    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss
     REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
     REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
     REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
-    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm
+    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
     REAL, DIMENSION(klon)              :: ycdragh, ycdragm
     REAL, DIMENSION(klon,klev)         :: yu, yv
@@ -470,8 +473,8 @@
 !****************************************************************************************
 ! Declarations specifiques pour le 1D. A reprendre 
+!****************************************************************************************
   REAL  :: fsens,flat
   LOGICAL :: ok_flux_surf ! initialized during first_call below
   COMMON /flux_arp/fsens,flat,ok_flux_surf
-!****************************************************************************************
 ! End of declarations
 !****************************************************************************************
@@ -546,11 +549,17 @@
     d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
     flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
-    d_u = 0.0     ; d_v = 0.0        ; yqsol = 0.0    
+    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0        ; yqsol = 0.0    
     ytherm = 0.0  ; ytke=0.
     
-    zcoefh(:,:) = 0.0
-    zcoefh(:,1) = 999999. ! zcoefh(:,k=1) should never be used
-    zcoefm(:,:) = 0.0
-    zcoefm(:,1) = 999999. !
+    tke(:,:,is_ave)=0.
+    IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
+       zcoefh(:,:,:) = 0.0
+       zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
+       zcoefm(:,:,:) = 0.0
+       zcoefm(:,1,:) = 999999. !
+    ELSE
+      zcoefm(:,:,is_ave)=0.
+      zcoefh(:,:,is_ave)=0.
+    ENDIF
     ytsoil = 999999. 
 
@@ -713,5 +722,5 @@
           ENDDO
        ENDDO
-       
+
        DO k = 1, nsoilmx
           DO j = 1, knon
@@ -747,4 +756,15 @@
             ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
             ycoefm, ycoefh, ytke)
+
+       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
+! In this case, coef_diff_turb is called for the Cd only
+       DO k = 2, klev
+          DO j = 1, knon
+             i = ni(j)
+             ycoefh(j,k)   = zcoefh(i,k,nsrf)
+             ycoefm(j,k)   = zcoefm(i,k,nsrf)
+          ENDDO
+       ENDDO
+       ENDIF
        
 !****************************************************************************************
@@ -924,4 +944,12 @@
 
 
+     y_d_t_diss(:,:)=0.
+     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
+        CALL yamada_c(knon,dtime,ypaprs,ypplay &
+    &   ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
+    &   ,iflag_pbl,nsrf)
+     ENDIF
+!     print*,'yamada_c OK'
+
        DO j = 1, knon
           y_dflux_t(j) = y_dflux_t(j) * ypct(j)
@@ -937,8 +965,8 @@
 !****************************************************************************************
 
-       tke(:,:,nsrf) = 0.
        DO k = 1, klev
           DO j = 1, knon
              i = ni(j)
+             y_d_t_diss(j,k)  = y_d_t_diss(j,k) * ypct(j)
              y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
              y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
@@ -951,8 +979,9 @@
              flux_v(i,k,nsrf) = y_flux_v(j,k)
 
-             tke(i,k,nsrf)    = ytke(j,k)
 
           ENDDO
        ENDDO
+
+!      print*,'Dans pbl OK1'
 
        evap(:,nsrf) = - flux_q(:,1,nsrf)
@@ -980,11 +1009,19 @@
        END DO
 
+!      print*,'Dans pbl OK2'
+
        DO k = 2, klev
           DO j = 1, knon
              i = ni(j)
-             zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)*ypct(j)
-             zcoefm(i,k) = zcoefm(i,k) + ycoefm(j,k)*ypct(j)
+             tke(i,k,nsrf)    = ytke(j,k)
+             zcoefh(i,k,nsrf) = ycoefh(j,k)
+             zcoefm(i,k,nsrf) = ycoefm(j,k)
+             tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
+             zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
+             zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
           END DO
        END DO
+
+!      print*,'Dans pbl OK3'
 
        IF ( nsrf .EQ. is_ter ) THEN 
@@ -1007,4 +1044,5 @@
           DO j = 1, knon
              i = ni(j)
+             d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
              d_t(i,k) = d_t(i,k) + y_d_t(j,k)
              d_q(i,k) = d_q(i,k) + y_d_q(j,k)
@@ -1013,4 +1051,6 @@
           END DO
        END DO
+
+!      print*,'Dans pbl OK4'
 
 !****************************************************************************************
@@ -1040,9 +1080,12 @@
 ! Calculations of diagnostic t,q at 2m and u, v at 10m
 
+!      print*,'Dans pbl OK41'
+!      print*,'tair1,yt(:,1),y_d_t(:,1)'
+!      print*, tair1,yt(:,1),y_d_t(:,1)
        DO j=1, knon
           i = ni(j)
           uzon(j) = yu(j,1) + y_d_u(j,1)
           vmer(j) = yv(j,1) + y_d_v(j,1)
-          tair1(j) = yt(j,1) + y_d_t(j,1)
+          tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
           qair1(j) = yq(j,1) + y_d_q(j,1)
           zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
@@ -1058,4 +1101,7 @@
        END DO
        
+!      print*,'Dans pbl OK42A'
+!      print*,'tair1,yt(:,1),y_d_t(:,1)'
+!      print*, tair1,yt(:,1),y_d_t(:,1)
 
 ! Calculate the temperature et relative humidity at 2m and the wind at 10m 
@@ -1064,4 +1110,5 @@
             tairsol, qairsol, rugo1, psfce, patm, &
             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+!      print*,'Dans pbl OK42B'
 
        DO j=1, knon
@@ -1077,4 +1124,5 @@
        END DO
 
+!      print*,'Dans pbl OK43'
 !IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
 !IM Ajoute dependance type surface
@@ -1093,4 +1141,5 @@
        END IF
 
+!   print*,'OK pbl 5'
        CALL HBTM(knon, ypaprs, ypplay, &
             yt2m,yt10m,yq2m,yq10m,yustar, &
@@ -1113,4 +1162,5 @@
        END DO
        
+!   print*,'OK pbl 6'
 #else 
 ! T2m not defined
@@ -1130,4 +1180,5 @@
 !****************************************************************************************
     
+!   print*,'OK pbl 7'
     zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
     zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
@@ -1143,4 +1194,5 @@
     END DO
 
+!   print*,'OK pbl 8'
     DO i = 1, klon
        zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
@@ -1161,4 +1213,5 @@
     s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
     
+!   print*,'OK pbl 9'
     
     DO nsrf = 1, nbsrf
@@ -1192,4 +1245,5 @@
        END DO
     END DO
+!   print*,'OK pbl 10'
 
     IF (check) THEN
Index: /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 1761)
@@ -55,6 +55,6 @@
       REAL, SAVE, ALLOCATABLE :: d_u_oli(:,:), d_v_oli(:,:)
       !$OMP THREADPRIVATE(d_u_oli, d_v_oli)
-      REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:)
-      !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf)
+      REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:), d_t_diss(:,:)
+      !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf,d_t_diss)
       REAL, SAVE, ALLOCATABLE :: d_u_vdf(:,:), d_v_vdf(:,:)
       !$OMP THREADPRIVATE(d_u_vdf, d_v_vdf)
@@ -216,5 +216,5 @@
       allocate(d_t_lscth(klon,klev),d_q_lscth(klon,klev))
       allocate(plul_st(klon),plul_th(klon))
-      allocate(d_t_vdf(klon,klev),d_q_vdf(klon,klev))
+      allocate(d_t_vdf(klon,klev),d_q_vdf(klon,klev),d_t_diss(klon,klev))
       allocate(d_u_vdf(klon,klev),d_v_vdf(klon,klev))
       allocate(d_t_oli(klon,klev),d_t_oro(klon,klev))
@@ -305,5 +305,5 @@
       deallocate(d_t_lscth,d_q_lscth)
       deallocate(plul_st,plul_th)
-      deallocate(d_t_vdf,d_q_vdf)
+      deallocate(d_t_vdf,d_q_vdf,d_t_diss)
       deallocate(d_u_vdf,d_v_vdf)
       deallocate(d_t_oli,d_t_oro)
Index: /LMDZ5/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1761)
@@ -156,4 +156,6 @@
   type(ctrl_out),save :: o_LWdnSFCclr   = ctrl_out((/ 1, 4, 10, 10, 5, 10 /),'LWdnSFCclr')
   type(ctrl_out),save :: o_bils         = ctrl_out((/ 1, 2, 10, 5, 10, 10 /),'bils')
+  type(ctrl_out),save :: o_bils_tke         = ctrl_out((/ 1, 2, 10, 5, 10, 10 /),'bils_tke')
+  type(ctrl_out),save :: o_bils_diss         = ctrl_out((/ 1, 2, 10, 5, 10, 10 /),'bils_diss')
   type(ctrl_out),save :: o_bils_ec         = ctrl_out((/ 1, 2, 10, 5, 10, 10 /),'bils_ec')
   type(ctrl_out),save :: o_bils_kinetic         = ctrl_out((/ 1, 2, 10, 5, 10, 10 /),'bils_kinetic')
@@ -554,4 +556,5 @@
   type(ctrl_out),save :: o_beta_prec    = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'beta_prec')
   type(ctrl_out),save :: o_dtvdf        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtvdf')
+  type(ctrl_out),save :: o_dtdis        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtdis')
   type(ctrl_out),save :: o_dqvdf        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqvdf')
   type(ctrl_out),save :: o_dteva        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dteva')
@@ -1043,4 +1046,6 @@
           CALL histdef2d(iff,clef_stations(iff),o_bils%flag,o_bils%name, "Surf. total heat flux", "W/m2")
           CALL histdef2d(iff,clef_stations(iff),o_bils_ec%flag,o_bils_ec%name, "Surf. total heat flux", "W/m2")
+          CALL histdef2d(iff,clef_stations(iff),o_bils_tke%flag,o_bils_tke%name, "Surf. total heat flux", "W/m2")
+          CALL histdef2d(iff,clef_stations(iff),o_bils_diss%flag,o_bils_diss%name, "Surf. total heat flux", "W/m2")
           CALL histdef2d(iff,clef_stations(iff),o_bils_kinetic%flag,o_bils_kinetic%name, "Surf. total heat flux", "W/m2")
           CALL histdef2d(iff,clef_stations(iff),o_bils_enthalp%flag,o_bils_enthalp%name, "Surf. total heat flux", "W/m2")
@@ -1625,4 +1630,5 @@
           CALL histdef3d(iff,clef_stations(iff),o_beta_prec%flag,o_beta_prec%name, "LS Conversion rate to prec", "(kg/kg)/s")
           CALL histdef3d(iff,clef_stations(iff),o_dtvdf%flag,o_dtvdf%name, "Boundary-layer dT", "K/s")
+          CALL histdef3d(iff,clef_stations(iff),o_dtdis%flag,o_dtdis%name, "TKE dissipation dT", "K/s")
           CALL histdef3d(iff,clef_stations(iff),o_dqvdf%flag,o_dqvdf%name, "Boundary-layer dQ", "(kg/kg)/s") 
           CALL histdef3d(iff,clef_stations(iff),o_dteva%flag,o_dteva%name, "Reevaporation dT", "K/s")
Index: /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 1761)
@@ -17,8 +17,11 @@
 !$OMP THREADPRIVATE(itau_con)
       REAL, ALLOCATABLE :: bils_ec(:) ! Contribution of energy conservation
+      REAL, ALLOCATABLE :: bils_tke(:) ! Contribution of energy conservation
+      REAL, ALLOCATABLE :: bils_diss(:) ! Contribution of energy conservation
       REAL, ALLOCATABLE :: bils_kinetic(:) ! bilan de chaleur au sol, kinetic
       REAL, ALLOCATABLE :: bils_enthalp(:) ! bilan de chaleur au sol
       REAL, ALLOCATABLE :: bils_latent(:) ! bilan de chaleur au sol
-!$OMP THREADPRIVATE(bils_ec,bils_kinetic,bils_enthalp,bils_latent)
+!$OMP THREADPRIVATE(bils_ec,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
+
 
 CONTAINS
@@ -32,5 +35,5 @@
       allocate(snow_o(klon), zfra_o(klon))
       allocate(itau_con(klon))
-      allocate (bils_ec(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
+      allocate (bils_ec(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
 
 END SUBROUTINE phys_output_var_init
@@ -42,5 +45,5 @@
 
       deallocate(snow_o,zfra_o,itau_con)
-      deallocate (bils_ec,bils_kinetic,bils_enthalp,bils_latent)
+      deallocate (bils_ec,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
 
 END SUBROUTINE phys_output_var_end
Index: /LMDZ5/trunk/libf/phylmd/phys_output_write.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1761)
@@ -360,7 +360,17 @@
        ENDIF
 
+       IF (o_bils_diss%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_bils_diss%name,itau_w,bils_diss)
+       ENDIF
+
        IF (o_bils_ec%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
      $o_bils_ec%name,itau_w,bils_ec)
+       ENDIF
+
+       IF (o_bils_tke%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_bils_tke%name,itau_w,bils_tke)
        ENDIF
 
@@ -1333,5 +1343,5 @@
 
 !====MS forcing diagnostics
-        if (new_aod) then	       
+        if (new_aod) then
         IF (o_swtoaas_nat%flag(iff)<=lev_files(iff)) THEN
         CALL histwrite_phy(nid_files(iff),clef_stations(iff),
@@ -1715,5 +1725,5 @@
        IF (o_kz%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_kz%name,itau_w,coefh)
+     $o_kz%name,itau_w,coefh(:,:,is_ave))
        ENDIF
 
@@ -1721,5 +1731,5 @@
        IF (o_kz_max%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_kz_max%name,itau_w,coefh)
+     $o_kz_max%name,itau_w,coefh(:,:,is_ave))
        ENDIF
       ENDIF
@@ -1908,4 +1918,10 @@
        ENDIF
 
+       IF (o_dtdis%flag(iff)<=lev_files(iff)) THEN
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_diss(1:klon,1:klev)/pdtphys
+      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
+     $o_dtdis%name,itau_w,zx_tmp_fi3d)
+       ENDIF
+
        IF (o_dqvdf%flag(iff)<=lev_files(iff)) THEN
       zx_tmp_fi3d(1:klon,1:klev)=d_q_vdf(1:klon,1:klev)/pdtphys
@@ -2191,5 +2207,5 @@
        IF (o_evu%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
-     $o_evu%name,itau_w,coefm)
+     $o_evu%name,itau_w,coefm(:,:,is_ave))
        ENDIF
 
Index: /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 1761)
@@ -63,5 +63,7 @@
 !$OMP THREADPRIVATE(ratqs)
       REAL, ALLOCATABLE, SAVE :: pbl_tke(:,:,:) ! turb kinetic energy
-!$OMP THREADPRIVATE(pbl_tke)
+      REAL, ALLOCATABLE, SAVE :: coefh(:,:,:) ! Kz enthalpie
+      REAL, ALLOCATABLE, SAVE :: coefm(:,:,:) ! Kz momentum
+!$OMP THREADPRIVATE(pbl_tke, coefh,coefm)
       REAL, ALLOCATABLE, SAVE :: zmax0(:), f0(:) ! 
 !$OMP THREADPRIVATE(zmax0,f0)
@@ -394,5 +396,7 @@
       ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
       ALLOCATE(ratqs(klon,klev))
-      ALLOCATE(pbl_tke(klon,klev+1,nbsrf))
+      ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1))
+      ALLOCATE(coefh(klon,klev+1,nbsrf+1))
+      ALLOCATE(coefm(klon,klev+1,nbsrf+1))
       ALLOCATE(zmax0(klon), f0(klon))
       ALLOCATE(ema_work1(klon,klev), ema_work2(klon,klev))
@@ -530,5 +534,5 @@
       deallocate(        u_ancien, v_ancien                 )
       deallocate(        tr_ancien)                           !RomP
-      deallocate(ratqs, pbl_tke)
+      deallocate(ratqs, pbl_tke,coefh,coefm)
       deallocate(zmax0, f0)
       deallocate(ema_work1, ema_work2)
Index: /LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1760)
+++ /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1761)
@@ -655,4 +655,5 @@
       real w0(klon)                                          ! Vitesse des thermiques au LCL
       real w_conv(klon)                                      ! Vitesse verticale de grande \'echelle au LCL
+      real tke0(klon,klev+1)                              ! TKE au début du pas de temps
       real therm_tke_max0(klon)                              ! TKE dans les thermiques au LCL 
       real env_tke_max0(klon)                                ! TKE dans l'environnement au LCL 
@@ -694,7 +695,4 @@
 cAA
 cAA  Pour phytrac 
-cAA
-      REAL coefh(klon,klev)     ! coef d'echange pour phytrac, valable pour 2<=k<=klev
-      REAL coefm(klon,klev)     ! coef d'echange pour U, V
       REAL u1(klon)             ! vents dans la premiere couche U
       REAL v1(klon)             ! vents dans la premiere couche V
@@ -1288,6 +1286,4 @@
           pbase=0
 cIM 180608
-c         pmflxr=0.
-c         pmflxs=0.
 
         itau_con=0
@@ -1396,4 +1392,9 @@
 
          CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
+         IF (klon_glo==1) THEN
+         coefh=0. ; coefm=0. ; pbl_tke=0.
+         coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
+         PRINT*,'FH WARNING : lignes a supprimer'
+         ENDIF
 cIM begin
           print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
@@ -1770,4 +1771,6 @@
       d1a(:,:)=0.
       dam(:,:)=0.
+          pmflxr=0.
+          pmflxs=0.
 ! RomP <<<
 
@@ -1785,4 +1788,5 @@
       ENDDO
       ENDDO
+      tke0(:,:)=pbl_tke(:,:,is_ave)
       IF (nqtot.GE.3) THEN
       DO iq = 3, nqtot
@@ -2073,5 +2077,5 @@
      s     albsol1,   albsol2,   sens,    evap,  
      s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
-     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
+     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss,
      s     coefh,     coefm,     slab_wfbils,                
      d     qsol,      zq2m,      s_pblh,  s_lcl,
@@ -2089,5 +2093,6 @@
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances de la diffusion turbulente
-      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
+      CALL add_phys_tend
+     s     (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
 !-----------------------------------------------------------------------------------------
 
@@ -3175,5 +3180,5 @@
      $                          paprs,
      $                          pplay,
-     $                          coefh,
+     $                          coefh(:,:,is_ave),
      $                          pphi,
      $                          t_seri,
@@ -3696,5 +3701,5 @@
      I     paprs,    pplay,     pmfu,     pmfd,
      I     pen_u,    pde_u,     pen_d,    pde_d,
-     I     cdragh,   coefh,     fm_therm, entr_therm,
+     I     cdragh,   coefh(:,:,is_ave),     fm_therm, entr_therm,
      I     u1,       v1,        ftsol,    pctsrf,
      I     ustar,     u10m,      v10m,
@@ -3723,5 +3728,5 @@
      I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      I                   fm_therm,entr_therm,
-     I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
+     I                   cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf,
      I                   frac_impa, frac_nucl,
      I                   pphis,airephy,dtime,itap,
@@ -3758,5 +3763,5 @@
       forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
       CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),
-     s        u_seri,v_seri,t_seri,q_seri,
+     s        u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:),
      s        zmasse,exner,d_t_ec)
       t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
@@ -3963,4 +3968,5 @@
 
 
+
 c=============================================================
 ! Separation entre thermiques et non thermiques dans les sorties
@@ -3991,5 +3997,4 @@
       endif
 
- 
 #include "phys_output_write.h"
 
@@ -4070,5 +4075,4 @@
          ENDIF !if callstats
      
-
       IF (lafin) THEN
          itau_phy = itau_phy + itap
Index: /LMDZ5/trunk/libf/phylmd/yamada_c.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/yamada_c.F90	(revision 1761)
+++ /LMDZ5/trunk/libf/phylmd/yamada_c.F90	(revision 1761)
@@ -0,0 +1,712 @@
+!
+! $Header$
+!
+      SUBROUTINE yamada_c(ngrid,timestep,plev,play &
+     &   ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar &
+     &   ,iflag_pbl,okiophys)
+      use dimphy
+      IMPLICIT NONE
+#include "iniprint.h"
+#include "YOMCST.h"
+!.......................................................................
+!ym#include "dimensions.h"
+!ym#include "dimphy.h"
+!.......................................................................
+!
+! timestep : pas de temps
+! g  : g
+! zlev : altitude a chaque niveau (interface inferieure de la couche
+!        de meme indice)
+! zlay : altitude au centre de chaque couche
+! u,v : vitesse au centre de chaque couche
+!       (en entree : la valeur au debut du pas de temps)
+! teta : temperature potentielle au centre de chaque couche
+!        (en entree : la valeur au debut du pas de temps)
+! cd : cdrag
+!      (en entree : la valeur au debut du pas de temps)
+! q2 : $q^2$ au bas de chaque couche
+!      (en entree : la valeur au debut du pas de temps)
+!      (en sortie : la valeur a la fin du pas de temps)
+! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
+!      couche)
+!      (en sortie : la valeur a la fin du pas de temps)
+! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
+!      (en sortie : la valeur a la fin du pas de temps)
+!
+!  iflag_pbl doit valoir entre 6 et 9
+!      l=6, on prend  systematiquement une longueur d'equilibre
+!    iflag_pbl=6 : MY 2.0
+!    iflag_pbl=7 : MY 2.0.Fournier
+!    iflag_pbl=8/9 : MY 2.5
+!       iflag_pbl=8 with special obsolete treatments for convergence
+!       with Cmpi5 NPv3.1 simulations
+!    iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
+!    iflag_pbl=12 = 11 with vertical diffusion off q2
+!
+!  2013/04/01 (FH hourdin@lmd.jussieu.fr)
+!     Correction for very stable PBLs (iflag_pbl=10 and 11)
+!     iflag_pbl=8 converges numerically with NPv3.1
+!     iflag_pbl=11 -> the model starts with NP from start files created by ce0l
+!                  -> the model can run with longer time-steps.
+!.......................................................................
+
+      REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t
+      REAL, DIMENSION(klon,klev) :: pu,pv,pt
+      REAL, DIMENSION(klon,klev) :: d_t_diss
+      INTEGER okiophys
+
+      REAL timestep
+      real plev(klon,klev+1)
+      real play(klon,klev)
+      real ustar(klon)
+      real kmin,qmin,pblhmin(klon),coriol(klon)
+      REAL zlev(klon,klev+1)
+      REAL zlay(klon,klev)
+      REAL zu(klon,klev)
+      REAL zv(klon,klev)
+      REAL zt(klon,klev)
+      REAL teta(klon,klev)
+      REAL cd(klon)
+      REAL q2(klon,klev+1),qpre
+      REAL unsdz(klon,klev)
+      REAL unsdzdec(klon,klev+1)
+
+      REAL km(klon,klev+1)
+      REAL kmpre(klon,klev+1),tmp2
+      REAL mpre(klon,klev+1)
+      REAL kn(klon,klev+1)
+      REAL kq(klon,klev+1)
+      real ff(klon,klev+1),delta(klon,klev+1)
+      real aa(klon,klev+1),aa0,aa1
+      integer iflag_pbl,ngrid
+      integer nlay,nlev
+
+      logical first
+      integer ipas
+      save first,ipas
+!FH/IM     data first,ipas/.true.,0/
+      data first,ipas/.false.,0/
+!$OMP THREADPRIVATE( first,ipas)
+
+      integer ig,k
+
+
+      real ri,zrif,zalpha,zsm,zsn
+      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
+
+      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
+      REAL, DIMENSION(klon,klev+1) :: km2,kn2,sqrtq
+      real dtetadz(klon,klev+1)
+      real m2cstat,mcstat,kmcstat
+      real l(klon,klev+1)
+      real leff(klon,klev+1)
+      real,allocatable,save :: l0(:)
+!$OMP THREADPRIVATE(l0)      
+      real sq(klon),sqz(klon),zz(klon,klev+1)
+      integer iter
+
+      real ric,rifc,b1,kap
+      save ric,rifc,b1,kap
+      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
+!$OMP THREADPRIVATE(ric,rifc,b1,kap)
+      real frif,falpha,fsm
+      real fl,zzz,zl0,zq2,zn2
+
+      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
+      real lyam(klon,klev),knyam(klon,klev)
+      real w2yam(klon,klev),t2yam(klon,klev)
+      logical,save :: firstcall=.true.
+
+REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
+REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
+REAL, DIMENSION(klon,klev) :: exner,masse
+REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
+
+!$OMP THREADPRIVATE(firstcall)       
+      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
+      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
+      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
+      fl(zzz,zl0,zq2,zn2)= &
+     &     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) &
+     &     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
+
+
+      if (firstcall) then
+        allocate(l0(klon))
+#ifdef IOPHYS
+        call iophys_ini
+#endif
+        firstcall=.false.
+      endif
+
+
+#ifdef IOPHYS
+if (okiophys==1) then
+call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
+call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
+endif
+#endif
+
+      nlay=klev
+      nlev=klev+1
+
+!-------------------------------------------------------------------------
+! Computation of conservative source terms from the turbulent tendencies
+!-------------------------------------------------------------------------
+
+
+   zu(:,:)=pu(:,:)+0.5*d_u(:,:)
+   zv(:,:)=pv(:,:)+0.5*d_v(:,:)
+   zt(:,:)=pt(:,:)+0.5*d_t(:,:)
+   do k=1,klev
+      exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA
+      masse(:,k)=(plev(:,k)-plev(:,k+1))/RG
+   enddo
+   teta(:,:)=zt(:,:)/exner(:,:)
+
+! Atmospheric mass at layer interfaces, where the TKE is computed
+   masseb(:,:)=0.
+   do k=1,klev
+      masseb(:,k)=masseb(:,k)+masse(:,k)
+      masseb(:,k+1)=masseb(:,k+1)+masse(:,k)
+    enddo
+    masseb(:,:)=0.5*masseb(:,:)
+
+
+
+   zlev(:,1)=0.
+   zlay(:,1)=RCPD*teta(:,1)*(1.-exner(:,1))
+   do k=1,klev-1
+      zlay(:,k+1)=zlay(:,k)+0.5*RCPD*(teta(:,k)+teta(:,k+1))*(exner(:,k)-exner(:,k+1))/RG
+      zlev(:,k)=0.5*(zlay(:,k)+zlay(:,k+1)) ! PASBO
+   enddo
+
+   fluxu(:,klev+1)=0.
+   fluxv(:,klev+1)=0.
+   fluxt(:,klev+1)=0.
+
+   do k=klev,1,-1
+      fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k)
+      fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k)
+      fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k) ! Flux de theta
+   enddo
+
+   dddu(:,1)=2*zu(:,1)*fluxu(:,1)
+   dddv(:,1)=2*zv(:,1)*fluxv(:,1)
+   dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1)
+
+   do k=2,klev
+      dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k)
+      dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k)
+      dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k)
+   enddo
+   dddu(:,klev+1)=0.
+   dddv(:,klev+1)=0.
+   dddt(:,klev+1)=0.
+
+#ifdef IOPHYS
+if (okiophys==1) then
+      call iophys_ecrit('zlay',klev,'Geop','m',zlay)
+      call iophys_ecrit('teta',klev,'teta','K',teta)
+      call iophys_ecrit('temp',klev,'temp','K',zt)
+      call iophys_ecrit('pt',klev,'temp','K',pt)
+      call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u)
+      call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v)
+      call iophys_ecrit('d_t',klev,'d_t','K/s',d_t)
+      call iophys_ecrit('exner',klev,'exner','',exner)
+      call iophys_ecrit('masse',klev,'masse','',masse)
+      call iophys_ecrit('masseb',klev,'masseb','',masseb)
+      call iophys_ecrit('Cm2',klev,'m2 conserv','m/s',(dddu(:,1:klev)+dddv(:,1:klev))/(masseb(:,1:klev)*timestep))
+      call iophys_ecrit('Cn2',klev,'m2 conserv','m/s',(rcpd*dddt(:,1:klev)/masseb(:,1:klev))/timestep)
+      call iophys_ecrit('rifc',klev,'rif conservative','',rcpd*dddt(:,1:klev)/min(dddu(:,1:klev)+dddv(:,1:klev),-1.e-20))
+endif
+#endif
+
+
+
+      ipas=ipas+1
+
+
+!.......................................................................
+!  les increments verticaux
+!.......................................................................
+!
+!!!!!! allerte !!!!!c
+!!!!!! zlev n'est pas declare a nlev !!!!!c
+!!!!!! ---->
+                                                      DO ig=1,ngrid
+            zlev(ig,nlev)=zlay(ig,nlay) &
+     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
+                                                      ENDDO
+!!!!!! <----
+!!!!!! allerte !!!!!c
+!
+      DO k=1,nlay
+                                                      DO ig=1,ngrid
+        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
+                                                      ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
+                                                      ENDDO
+      DO k=2,nlay
+                                                      DO ig=1,ngrid
+        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
+                                                     ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
+                                                     ENDDO
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Computing M^2, N^2, Richardson numbers, stability functions
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+      do k=2,klev
+                                                          do ig=1,ngrid
+         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
+         m2(ig,k)=((zu(ig,k)-zu(ig,k-1))**2+(zv(ig,k)-zv(ig,k-1))**2)/(dz(ig,k)*dz(ig,k))
+         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
+         n2(ig,k)=RG*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
+!        n2(ig,k)=0.
+         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
+         if (ri.lt.ric) then
+            rif(ig,k)=frif(ri)
+         else
+            rif(ig,k)=rifc
+         endif
+         if(rif(ig,k).lt.0.16) then
+            alpha(ig,k)=falpha(rif(ig,k))
+            sm(ig,k)=fsm(rif(ig,k))
+         else
+            alpha(ig,k)=1.12
+            sm(ig,k)=0.085
+         endif
+         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
+                                                          enddo
+      enddo
+
+
+
+!====================================================================
+!  Computing the mixing length
+!====================================================================
+
+!   Mise a jour de l0
+      if (iflag_pbl==8.or.iflag_pbl==10) then
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Iterative computation of l0
+! This version is kept for iflag_pbl only for convergence
+! with NPv3.1 Cmip5 simulations
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+                                                          do ig=1,ngrid
+      sq(ig)=1.e-10
+      sqz(ig)=1.e-10
+                                                          enddo
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        zq=sqrt(q2(ig,k))
+        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+      enddo
+                                                          do ig=1,ngrid
+      l0(ig)=0.2*sqz(ig)/sq(ig)
+                                                          enddo
+      do k=2,klev
+                                                          do ig=1,ngrid
+         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+                                                          enddo
+      enddo
+!     print*,'L0 cas 8 ou 10 ',l0
+
+      else
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! In all other case, the assymptotic mixing length l0 is imposed (100m)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+          l0(:)=150.
+          do k=2,klev
+                                                          do ig=1,ngrid
+             l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+                                                          enddo
+          enddo
+!     print*,'L0 cas autres ',l0
+
+      endif
+
+
+#ifdef IOPHYS
+if (okiophys==1) then
+call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
+call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
+call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
+call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
+endif
+#endif
+
+
+IF (iflag_pbl<20) then
+      ! For diagnostics only
+      RETURN
+
+ELSE
+
+!  print*,'OK1'
+
+! Evolution of TKE under source terms K M2 and K N2
+   leff(:,:)=max(l(:,:),1.)
+   IF (iflag_pbl==29) THEN
+      km2(:,:)=km(:,:)*m2(:,:)
+      kn2(:,:)=kn2(:,:)*rif(:,:)
+   ELSEIF (iflag_pbl==25) THEN
+      DO k=1,klev
+         km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
+         &        /(masse(:,k)*timestep)
+         kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
+         leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
+      ENDDO
+      km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
+   ELSE
+      km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep)
+      kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep)
+   ENDIF
+   q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:))
+   q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4)
+
+! Dissipation of TKE
+   q2old(:,:)=q2(:,:)
+   q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1))
+   q2(:,:)=q2(:,:)*q2(:,:)
+   IF (iflag_pbl<=24) THEN
+      DO k=1,klev
+         d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k))
+      ENDDO
+   ELSE IF (iflag_pbl<=27) THEN
+      DO k=1,klev
+         d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
+      ENDDO
+   ENDIF
+!  print*,'iflag_pbl ',d_t_diss
+
+
+! Compuation of stability functions
+   IF (iflag_pbl/=29) THEN
+      DO k=1,klev
+      DO ig=1,ngrid
+         IF (ABS(km2(ig,k))<=1.e-20) THEN
+            rif(ig,k)=0.
+         ELSE
+            rif(ig,k)=min(kn2(ig,k)/km2(ig,k),rifc)
+         ENDIF
+         IF (rif(ig,k).lt.0.16) THEN
+            alpha(ig,k)=falpha(rif(ig,k))
+            sm(ig,k)=fsm(rif(ig,k))
+         else
+            alpha(ig,k)=1.12
+            sm(ig,k)=0.085
+         endif
+      ENDDO
+      ENDDO
+    ENDIF
+
+! Computation of turbulent diffusivities
+   IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
+     DO k=2,klev
+        sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
+     ENDDO
+   ELSE
+     DO k=2,klev
+        sqrtq(:,k)=sqrt(q2(:,k))
+     ENDDO
+   ENDIF
+   DO k=2,klev
+   DO ig=1,ngrid
+      km(ig,k)=leff(ig,k)*sqrtq(ig,k)*sm(ig,k)
+      kn(ig,k)=km(ig,k)*alpha(ig,k)
+      kq(ig,k)=leff(ig,k)*zq*0.2
+!         print*,q2(ig,k),zq,km(ig,k)
+   ENDDO
+   ENDDO
+
+
+
+#ifdef IOPHYS
+if (okiophys==1) then
+call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
+call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
+call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
+call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
+call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
+call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
+call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
+call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
+call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
+call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
+endif
+#endif
+
+ENDIF
+
+
+!  print*,'OK2'
+      RETURN
+!====================================================================
+!   Yamada 2.0
+!====================================================================
+      if (iflag_pbl.eq.6) then
+
+      do k=2,klev
+         q2(:,k)=l(:,k)**2*zz(:,k)
+      enddo
+
+
+      else if (iflag_pbl.eq.7) then
+!====================================================================
+!   Yamada 2.Fournier
+!====================================================================
+
+!  Calcul de l,  km, au pas precedent
+      do k=2,klev
+                                                          do ig=1,ngrid
+!        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         mpre(ig,k)=sqrt(m2(ig,k))
+!        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+                                                          enddo
+      enddo
+
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
+        mcstat=sqrt(m2cstat)
+
+!        print*,'M2 L=',k,mpre(ig,k),mcstat
+!
+!  -----{puis on ecrit la valeur de q qui annule l'equation de m
+!        supposee en q3}
+!
+        IF (k.eq.2) THEN
+          kmcstat=1.E+0 / mcstat &
+     &    *( unsdz(ig,k)*kmpre(ig,k+1) &
+     &                        *mpre(ig,k+1) &
+     &      +unsdz(ig,k-1) &
+     &              *cd(ig) &
+     &              *( sqrt(zu(ig,3)**2+zv(ig,3)**2) &
+     &                -mcstat/unsdzdec(ig,k) &
+     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) &
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ELSE
+          kmcstat=1.E+0 / mcstat &
+     &    *( unsdz(ig,k)*kmpre(ig,k+1) &
+     &                        *mpre(ig,k+1) &
+     &      +unsdz(ig,k-1)*kmpre(ig,k-1) &
+     &                          *mpre(ig,k-1) ) &
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ENDIF
+!       print*,'T2 L=',k,tmp2
+        tmp2=kmcstat &
+     &      /( sm(ig,k)/q2(ig,k) ) &
+     &      /l(ig,k)
+        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
+!       print*,'Q2 L=',k,q2(ig,k)
+!
+                                                          enddo
+      enddo
+
+      else if (iflag_pbl==8.or.iflag_pbl==9) then
+!====================================================================
+!   Yamada 2.5 a la Didi
+!====================================================================
+
+
+!  Calcul de l,  km, au pas precedent
+      do k=2,klev
+                                                          do ig=1,ngrid
+!        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         if (delta(ig,k).lt.1.e-20) then
+!     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
+            delta(ig,k)=1.e-20
+         endif
+         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         aa0=(m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
+         aa1=(m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
+! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
+         aa(ig,k)=aa1*timestep/(delta(ig,k)*l(ig,k))
+!     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+         qpre=sqrt(q2(ig,k))
+!        if (iflag_pbl.eq.8 ) then
+            if (aa(ig,k).gt.0.) then
+               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
+            else
+               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+            endif
+!        else ! iflag_pbl=9
+!           if (aa(ig,k)*qpre.gt.0.9) then
+!              q2(ig,k)=(qpre*10.)**2
+!           else
+!              q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+!           endif
+!        endif
+         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
+!     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
+                                                          enddo
+      enddo
+
+      else if (iflag_pbl>=10) then
+
+!        print*,'Schema mixte D'
+!        print*,'Longueur ',l(:,:)
+         do k=2,klev-1
+            l(:,k)=max(l(:,k),1.)
+            km(:,k)=l(:,k)*sqrt(q2(:,k))*sm(:,k)
+            q2(:,k)=q2(:,k)+timestep*km(:,k)*m2(:,k)*(1.-rif(:,k))
+            q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
+            q2(:,k)=1./(1./sqrt(q2(:,k))+timestep/(2*l(:,k)*b1))
+            q2(:,k)=q2(:,k)*q2(:,k)
+         enddo
+
+
+      else
+         stop'Cas nom prevu dans yamada4'
+
+      endif ! Fin du cas 8
+
+!     print*,'OK8'
+
+!====================================================================
+!   Calcul des coefficients de m�ange
+!====================================================================
+      do k=2,klev
+!     print*,'k=',k
+                                                          do ig=1,ngrid
+!abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
+         zq=sqrt(q2(ig,k))
+         km(ig,k)=l(ig,k)*zq*sm(ig,k)
+         kn(ig,k)=km(ig,k)*alpha(ig,k)
+         kq(ig,k)=l(ig,k)*zq*0.2
+!     print*,'KML=',km(ig,k),kn(ig,k)
+                                                          enddo
+      enddo
+
+! Transport diffusif vertical de la TKE.
+      if (iflag_pbl.ge.12) then
+!       print*,'YAMADA VDIF'
+        q2(:,1)=q2(:,2)
+        call vdif_q2(timestep,RG,RD,ngrid,plev,zt,kq,q2)
+      endif
+
+!   Traitement des cas noctrunes avec l'introduction d'une longueur
+!   minilale.
+
+!====================================================================
+!   Traitement particulier pour les cas tres stables.
+!   D'apres Holtslag Boville.
+
+      if (prt_level>1) THEN
+       print*,'YAMADA4 0'
+      endif !(prt_level>1) THEN
+                                                          do ig=1,ngrid
+      coriol(ig)=1.e-4
+      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
+                                                          enddo
+
+!     print*,'pblhmin ',pblhmin
+!Test a remettre 21 11 02
+! test abd 13 05 02      if(0.eq.1) then
+      if(1==1) then
+      if(iflag_pbl==8.or.iflag_pbl==10) then
+
+      do k=2,klev
+         do ig=1,ngrid
+            if (teta(ig,2).gt.teta(ig,1)) then
+               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
+               kmin=kap*zlev(ig,k)*qmin
+            else
+               kmin=-1. ! kmin n'est utilise que pour les SL stables.
+            endif 
+            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
+!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
+!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
+               kn(ig,k)=kmin
+               km(ig,k)=kmin
+               kq(ig,k)=kmin
+!   la longueur de melange est suposee etre l= kap z
+!   K=l q Sm d'ou q2=(K/l Sm)**2
+               q2(ig,k)=(qmin/sm(ig,k))**2
+            endif
+         enddo
+      enddo
+
+      else
+
+      do k=2,klev
+         do ig=1,ngrid
+            if (teta(ig,2).gt.teta(ig,1)) then
+               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
+               kmin=kap*zlev(ig,k)*qmin
+            else
+               kmin=-1. ! kmin n'est utilise que pour les SL stables.
+            endif 
+            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
+!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
+!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
+               kn(ig,k)=kmin
+               km(ig,k)=kmin
+               kq(ig,k)=kmin
+!   la longueur de melange est suposee etre l= kap z
+!   K=l q Sm d'ou q2=(K/l Sm)**2
+               sm(ig,k)=1.
+               alpha(ig,k)=1.
+               q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
+               zq=sqrt(q2(ig,k))
+               km(ig,k)=l(ig,k)*zq*sm(ig,k)
+               kn(ig,k)=km(ig,k)*alpha(ig,k)
+               kq(ig,k)=l(ig,k)*zq*0.2
+            endif
+         enddo
+      enddo
+      endif
+
+      endif
+
+      if (prt_level>1) THEN
+       print*,'YAMADA4 1'
+      endif !(prt_level>1) THEN
+!   Diagnostique pour stokage
+
+      if(1.eq.0)then
+      rino=rif
+      smyam(1:ngrid,1)=0.
+      styam(1:ngrid,1)=0.
+      lyam(1:ngrid,1)=0.
+      knyam(1:ngrid,1)=0.
+      w2yam(1:ngrid,1)=0.
+      t2yam(1:ngrid,1)=0.
+
+      smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
+      styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
+      lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
+      knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
+
+!   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
+
+      w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 &
+     &    +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev) &
+     &    *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
+
+      t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) &
+     &    *dtetadz(1:ngrid,2:klev)**2 &
+     &    /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
+      endif
+
+!     print*,'OKFIN'
+      first=.false.
+      return
+      end
