Index: /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4881)
@@ -14,5 +14,5 @@
   SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
        ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
-       ycoefm, ycoefh ,yq2, ydrgpro)
+       ycoefm, ycoefh ,yq2, yeps, ydrgpro)
  
     USE dimphy
@@ -47,4 +47,5 @@
 ! Output arguments
 !****************************************************************************************
+    REAL, DIMENSION(klon,klev+1), INTENT(OUT)  :: yeps
     REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefh
     REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
@@ -69,5 +70,6 @@
     ykmq = 0 !ym missing init
     
-    
+    yeps(:,:) = 0.
+
 !****************************************************************************************    
 ! Calcul de coefficients de diffusion turbulent de l'atmosphere : 
@@ -182,5 +184,5 @@
           CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &
                yzlev,yzlay,yu,yv,yteta, &
-               ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &
+               ycdragm,yq2,yeps,ykmm,ykmn,ykmq,yustar, &
                iflag_pbl,ydrgpro)
 !FC
Index: /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4881)
@@ -7,5 +7,5 @@
 subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
                         wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
-                        tke,Km_out,Kh_out)
+                        tke,eps,Km_out,Kh_out)
 
 !========================================================================
@@ -70,14 +70,15 @@
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
-REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
-REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
-REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
-REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
-REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv   ! surface drag coefficient for momentum
-
-REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
-
-REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
-REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp     ! temperature (K)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap     ! specific humidity (kg/kg)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play     ! pressure (Pa)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf  ! pressure at interfaces(Pa)
+REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv ! surface drag coefficient for momentum
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke      ! turbulent kinetic energy at interface between layers (m2/s2)
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: eps      ! output: TKE dissipation rate at interface between layers (m2/s3)
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers (m2/s)
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers (m2/s)
 
 ! Local variables
@@ -259,4 +260,5 @@
                 tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
                 shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
+                eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
             ENDDO
         ENDDO
@@ -275,5 +277,6 @@
             qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
             qq=max(0.,qq)
-            tke(igrid,ilay)=0.5*(qq**2) 
+            tke(igrid,ilay)=0.5*(qq**2)
+            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) 
             ENDDO
         ENDDO
@@ -292,4 +295,5 @@
                 /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
             tke(igrid,ilay)=0.5*(qq**2) 
+            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
             ENDDO
         ENDDO
@@ -302,4 +306,5 @@
         DO ilay=2,nlay
             DO igrid=1,ngrid
+            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
             qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
             IF (Ri(igrid,ilay) .LT. 0.) THEN
@@ -327,4 +332,5 @@
             qq=max(0.,qq)
             tke(igrid,ilay)=0.5*(qq**2)
+            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
             ENDDO
         ENDDO
Index: /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 4881)
@@ -10,5 +10,5 @@
 subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
                         wind_u,wind_v,temp,qvap,play,pinterf, &
-                        tke,Km_out,Kh_out)
+                        tke,eps,Km_out,Kh_out)
 
 
@@ -43,6 +43,7 @@
 
 
-REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke      ! turbulent kinetic energy at interface between layers
 
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: eps      ! output: tke dissipation rate at interface between layers
 REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
 REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
@@ -56,5 +57,5 @@
 call atke_compute_km_kh(ngrid,nlay,dtime,&
                         wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
-                        tke,Km_out,Kh_out)
+                        tke,eps,Km_out,Kh_out)
 
 
@@ -75,5 +76,5 @@
    call atke_compute_km_kh(ngrid,nlay,dtime,&
                         wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
-                        tke,Km_out,Kh_out)
+                        tke,eps,Km_out,Kh_out)
 
 end if
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4881)
@@ -231,5 +231,5 @@
 !jyg<
 !!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
-       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x,  &
+       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x, eps_x, &
 !>jyg
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -408,4 +408,6 @@
 ! Output variables
 !****************************************************************************************
+    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(OUT)   :: eps_x      ! TKE dissipation rate
+
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
@@ -627,5 +629,5 @@
     REAL, DIMENSION(klon,klev)         :: delp
     REAL, DIMENSION(klon,klev+1)       :: ypaprs
-    REAL, DIMENSION(klon,klev+1)       :: ytke
+    REAL, DIMENSION(klon,klev+1)       :: ytke, yeps
     REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
 !FC 
@@ -647,5 +649,5 @@
     REAL, DIMENSION(klon)              :: ywake_s, ywake_cstar, ywake_dens
 !!!
-    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w
+    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w, yeps_x, yeps_w
     REAL, DIMENSION(klon,klev+1)       :: ywake_dltke
     REAL, DIMENSION(klon,klev)         :: yu_x, yv_x, yu_w, yv_w
@@ -1046,4 +1048,5 @@
 !!    tke(:,:,is_ave)=0.
     tke_x(:,:,is_ave)=0.
+    eps_x(:,:,is_ave)=0.
 
     wake_dltke(:,:,is_ave)=0.
@@ -1080,4 +1083,5 @@
 
     ytke=0.
+    yeps=0.
     yri0(:)=0.
 !FC 
@@ -1093,4 +1097,5 @@
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
     ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
+    yeps_x=0.     ; yeps_w=0.
     y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
 !!    d_t_w=0.      ; d_q_w=0.         
@@ -1701,5 +1706,5 @@
         CALL call_atke(dtime,knon,klev,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
                   yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),       &
-                  ytke(1:knon,:),ycoefm(1:knon,:), ycoefh(1:knon,:))
+                  ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
 
         ELSE
@@ -1707,5 +1712,5 @@
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
-            ycoefm, ycoefh, ytke, y_treedrg)
+            ycoefm, ycoefh, ytke, yeps, y_treedrg)
 !            ycoefm, ycoefh, ytke)
 !FC y_treedrg ajoute
@@ -1746,5 +1751,5 @@
         CALL call_atke(dtime,knon,klev,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
                        yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),  &
-                       ytke_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
+                       ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
 
         ELSE
@@ -1752,5 +1757,5 @@
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
-            ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
+            ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)
 !            ycoefm_x, ycoefh_x, ytke_x)
 !FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
@@ -1786,5 +1791,5 @@
         CALL call_atke(dtime,knon,klev,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
                 yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),      &
-                ytke_w(1:knon,:),ycoefm_w(1:knon,:), ycoefh_w(1:knon,:))
+                ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
 
         ELSE
@@ -1792,5 +1797,5 @@
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
-            ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
+            ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)
 !            ycoefm_w, ycoefh_w, ytke_w)
        IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
@@ -2855,6 +2860,7 @@
 !!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
               tke_x(i,k,nsrf)    = ytke(j,k)
-              tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
-
+              tke_x(i,k,is_ave)  = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
+              eps_x(i,k,nsrf)    = yeps(j,k)
+              eps_x(i,k,is_ave)  = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)
 !>jyg
            ENDDO
@@ -2870,5 +2876,7 @@
 !!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
             tke_x(i,k,nsrf)   = ytke_x(j,k)
-            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)        
+            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)       
+            eps_x(i,k,nsrf)   = yeps_x(j,k)
+            eps_x(i,k,is_ave)   = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j) 
             wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
             
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4881)
@@ -22,6 +22,8 @@
       REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
       !$OMP THREADPRIVATE(d_rneb_dyn)
-      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),tke_dissip(:,:,:),wprime(:,:,:)
-      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip,wprime)
+      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
+      !$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
+      REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
+      !$OMP THREADPRIVATE(pbl_eps)
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -666,7 +668,7 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
-      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),tke_dissip(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
-      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;tke_dissip(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
-
+      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
+      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
       ALLOCATE(rhcl(klon,klev))
       ALLOCATE(tr_seri(klon,klev,nbtr))
@@ -1019,6 +1021,6 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
       DEALLOCATE(u_seri,v_seri)
-      DEALLOCATE(l_mixmin,l_mix, tke_dissip,wprime)
-
+      DEALLOCATE(l_mixmin,l_mix,wprime)
+      DEALLOCATE(pbl_eps)
       DEALLOCATE(rhcl)
       DEALLOCATE(tr_seri)
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4881)
@@ -311,5 +311,5 @@
          zn2mout, t2m_min_mon, t2m_max_mon, evap, &
          snowerosion, zxustartlic, zxrhoslic, zxqsaltlic, &
-         l_mixmin,l_mix, tke_dissip, &
+         l_mixmin,l_mix, pbl_eps, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
          rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
@@ -1302,5 +1302,5 @@
                 DO k=1,klev
                    zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
-                        +pctsrf(:,nsrf)*tke_dissip(:,k,nsrf)
+                        +pctsrf(:,nsrf)*pbl_eps(:,k,nsrf)
                 ENDDO
              ENDDO
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4881)
@@ -282,5 +282,5 @@
        zxfluxlat_x, zxfluxlat_w, &
        !
-       pbl_tke_input, tke_dissip, l_mix, wprime,&
+       pbl_tke_input, pbl_eps, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
        cdragh_x, cdragh_w, &
@@ -1257,5 +1257,5 @@
     integer iostat
 
-    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
+    REAL, dimension(klon,klev+1) :: l_mix_ave, wprime_ave
     REAL zzz
     !albedo SB >>>
@@ -2905,5 +2905,5 @@
             fluxt,   fluxu,  fluxv, &
             dsens,     devap,     zxsnow, &
-            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, &
+            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, pbl_eps,  &
                                 !nrlmd+jyg<
             wake_delta_pbl_TKE, &
@@ -3849,5 +3849,4 @@
     ! water distribution
 
-    tke_dissip_ave(:,:)=0.
     l_mix_ave(:,:)=0.
     wprime_ave(:,:)=0.
@@ -3855,5 +3854,4 @@
     DO nsrf = 1, nbsrf
        DO i = 1, klon
-          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
           l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
           wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
@@ -3869,5 +3867,5 @@
          paprs,pplay,t_seri,q_seri, &
          qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
-         omega,pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave, &
+         omega,pbl_tke(:,:,is_ave),pbl_eps(:,:,is_ave),l_mix_ave,wprime_ave, &
          t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
          ratqs,ratqsc,ratqs_inter_)
Index: /LMDZ6/trunk/libf/phylmd/yamada4.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 4881)
@@ -2,8 +2,8 @@
 
 SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
-    cd, tke, km, kn, kq, ustar, iflag_pbl, drgpro)
+    cd, tke, eps, km, kn, kq, ustar, iflag_pbl, drgpro)
 
   USE dimphy, only : klev,klon
-  USE phys_local_var_mod, only: tke_dissip,wprime
+  USE phys_local_var_mod, only: wprime
   USE yamada_ini_mod, only : new_yamada4,yamada4_num,hboville
   USE yamada_ini_mod, only : prt_level, lunout,pbl_lmixmin_alpha,b1,kap,viscom,viscoh
@@ -69,4 +69,5 @@
   ! Outputs:
   !==========
+  ! eps: tke dissipation rate 
   ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
   ! couche)
@@ -95,4 +96,5 @@
   REAL cd(klon)
   REAL tke(klon, klev+1)
+  REAL eps(klon,klev+1)
   REAL unsdz(klon, klev)
   REAL unsdzdec(klon, klev+1)
@@ -686,5 +688,6 @@
 
 ! Diagnostics
- tke_dissip(1:ngrid,:,nsrf)=0.
+ 
+ eps(:,:)=dissip(:,:)
  wprime(1:ngrid,:,nsrf)=0.
  DO k=2,klev
@@ -692,7 +695,7 @@
        jg=ni(ig)
        wprime(jg,k,nsrf)=sqrt(MAX(1./3*q2(ig,k),0.))
-       tke_dissip(jg,k,nsrf)=dissip(ig,k)
     ENDDO
  ENDDO
+
  
 !=============================================================================
Index: /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4881)
@@ -24,5 +24,5 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
-  USE lmdz_call_atke,      ONLY :  call_atke
+  USE lmdz_call_atke,      ONLY : call_atke
   USE ioipsl_getin_p_mod,  ONLY : getin_p
   USE cdrag_mod
@@ -308,5 +308,5 @@
 !jyg<
 !!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
-       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x,  &
+       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x, eps_x, &
 !>jyg
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -506,4 +506,6 @@
 ! Output variables
 !****************************************************************************************
+    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(OUT)   :: eps_x      ! TKE dissipation rate
+
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
@@ -765,5 +767,5 @@
     REAL, DIMENSION(klon,klev)         :: delp
     REAL, DIMENSION(klon,klev+1)       :: ypaprs
-    REAL, DIMENSION(klon,klev+1)       :: ytke
+    REAL, DIMENSION(klon,klev+1)       :: ytke, yeps
     REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
 !FC 
@@ -785,5 +787,5 @@
     REAL, DIMENSION(klon)              :: ywake_s, ywake_cstar, ywake_dens
 !!!
-    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w
+    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w, yeps_x, yeps_w
     REAL, DIMENSION(klon,klev+1)       :: ywake_dltke
     REAL, DIMENSION(klon,klev)         :: yu_x, yv_x, yu_w, yv_w
@@ -1261,4 +1263,5 @@
 !!    tke(:,:,is_ave)=0.
     tke_x(:,:,is_ave)=0.
+    eps_x(:,:,is_ave)=0.
 
     wake_dltke(:,:,is_ave)=0.
@@ -1295,4 +1298,5 @@
 
     ytke=0.
+    yeps=0.
     yri0(:)=0.
 !FC 
@@ -1308,4 +1312,5 @@
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
     ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
+    yeps_x=0.     ; yeps_w=0.
     y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
 !!    d_t_w=0.      ; d_q_w=0.         
@@ -1979,5 +1984,5 @@
         CALL call_atke(dtime,knon,klev,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
                   yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),       &
-                  ytke(1:knon,:),ycoefm(1:knon,:), ycoefh(1:knon,:))
+                  ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
 
         ELSE
@@ -1985,5 +1990,5 @@
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
-            ycoefm, ycoefh, ytke, y_treedrg)
+            ycoefm, ycoefh, ytke, yeps, y_treedrg)
 !            ycoefm, ycoefh, ytke)
 !FC y_treedrg ajoute
@@ -2023,5 +2028,5 @@
         CALL call_atke(dtime,knon,klev,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
                        yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),  &
-                       ytke_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
+                       ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
 
         ELSE
@@ -2029,5 +2034,5 @@
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
-            ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
+            ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)
 !            ycoefm_x, ycoefh_x, ytke_x)
 !FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
@@ -2063,12 +2068,11 @@
         CALL call_atke(dtime,knon,klev,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
                 yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),      &
-                ytke_w(1:knon,:),ycoefm_w(1:knon,:), ycoefh_w(1:knon,:))
+                ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
 
         ELSE
-
 
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
-            ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
+            ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)
 !            ycoefm_w, ycoefh_w, ytke_w)
        IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
@@ -3332,6 +3336,7 @@
 !!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
               tke_x(i,k,nsrf)    = ytke(j,k)
-              tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
-
+              tke_x(i,k,is_ave)  = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
+              eps_x(i,k,nsrf)    = yeps(j,k)
+              eps_x(i,k,is_ave)  = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)
 !>jyg
            ENDDO
@@ -3347,5 +3352,7 @@
 !!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
             tke_x(i,k,nsrf)   = ytke_x(j,k)
-            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)        
+            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)       
+            eps_x(i,k,nsrf)   = yeps_x(j,k)
+            eps_x(i,k,is_ave)   = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j) 
             wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
             
Index: /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4881)
@@ -22,6 +22,8 @@
       REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
       !$OMP THREADPRIVATE(d_rneb_dyn)
-      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),tke_dissip(:,:,:),wprime(:,:,:)
-      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip,wprime)
+      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
+      !$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
+      REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
+      !$OMP THREADPRIVATE(pbl_eps)
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -764,7 +766,7 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
-      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),tke_dissip(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
-      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;tke_dissip(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
-
+      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
+      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
+      ALLOCATE(pbl_eps(klon,klev+1,nbsrf))      
       ALLOCATE(tr_seri(klon,klev,nbtr))
       ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
@@ -1175,5 +1177,6 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
       DEALLOCATE(u_seri,v_seri)
-      DEALLOCATE(l_mixmin,l_mix, tke_dissip,wprime)
+      DEALLOCATE(l_mixmin,l_mix,wprime)
+      DEALLOCATE(pbl_eps)
 
       DEALLOCATE(tr_seri)
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4880)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4881)
@@ -294,5 +294,5 @@
        zxfluxlat_x, zxfluxlat_w, &
        !
-       pbl_tke_input, tke_dissip, l_mix, wprime,&
+       pbl_tke_input, pbl_eps, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
        cdragh_x, cdragh_w, &
@@ -1321,5 +1321,5 @@
     integer iostat
 
-    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
+    REAL, dimension(klon,klev+1) :: l_mix_ave, wprime_ave
     REAL zzz
     !albedo SB >>>
@@ -3192,5 +3192,5 @@
             fluxt,   fluxu,  fluxv, &
             dsens,     devap,     zxsnow, &
-            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, &
+            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, pbl_eps,  &
                                 !nrlmd+jyg<
             wake_delta_pbl_TKE, &
@@ -4756,5 +4756,4 @@
     ! water distribution
 
-    tke_dissip_ave(:,:)=0.
     l_mix_ave(:,:)=0.
     wprime_ave(:,:)=0.
@@ -4762,5 +4761,4 @@
     DO nsrf = 1, nbsrf
        DO i = 1, klon
-          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
           l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
           wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
@@ -4776,5 +4774,5 @@
          paprs,pplay,t_seri,q_seri, &
          qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
-         omega,pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave, &
+         omega,pbl_tke(:,:,is_ave),pbl_eps(:,:,is_ave),l_mix_ave,wprime_ave, &
          t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
          ratqs,ratqsc,ratqs_inter_)
