Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 659)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 660)
@@ -6,5 +6,6 @@
      & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
      & pplay,pplev,pphi,zpopsk, &
-     & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar)
+     & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke, &
+     & pdhdif,hfmax,wstar,sensibFlux)
 
        USE ioipsl_getincom
@@ -34,4 +35,6 @@
       REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
       REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pdhdif(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: sensibFlux(ngridmx)
 
 !--------------------------------------------------------
@@ -66,5 +69,5 @@
       REAL masse(ngridmx,nlayermx)
       LOGICAL qtransport_thermals,dtke_thermals
-      INTEGER l,ig,iq,ii(1)
+      INTEGER l,ig,iq,ii(1),k
       CHARACTER (LEN=20) modname
 
@@ -100,4 +103,7 @@
       REAL hfmax(ngridmx),wmax(ngridmx)
       REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
+      REAL rpdhd(ngridmx,nlayermx)
+      REAL wtdif(ngridmx,nlayermx),rho(ngridmx,nlayermx)
+      REAL wtth(ngridmx,nlayermx)
 
 !--------------------------------------------------------
@@ -380,5 +386,4 @@
 
       DO ig=1,ngridmx
-         hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
          wmax(ig)=MAXVAL(zw2(ig,:))
       ENDDO
@@ -454,4 +459,57 @@
       ENDDO
 
+! we must add the heat flux from the diffusion scheme to hfmax
+
+! compute rho as it is after the diffusion
+
+      rho(:,:)=pplay(:,:)                                               &
+     & /(r*(pt(:,:)+pdhdif(:,:)*zpopsk(:,:)*ptimestep))
+
+! integrate -rho*pdhdif
+
+      rpdhd(:,:)=0.
+
+      DO ig=1,ngridmx
+       DO l=1,lmax(ig)
+        rpdhd(ig,l)=0.
+        DO k=1,l
+         rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*pdhdif(ig,k)*                &
+     & (zzlev(ig,k+1)-zzlev(ig,k))
+        ENDDO
+        rpdhd(ig,l)=rpdhd(ig,l)-sensibFlux(ig)/cpp
+       ENDDO
+      ENDDO
+
+! compute w'teta' from diffusion
+
+      wtdif(:,:)=rpdhd(:,:)/rho(:,:)
+
+! compute rho as it is after the thermals
+
+      rho(:,:)=pplay(:,:)/(r*(zt(:,:)))
+! integrate -rho*pdhdif
+
+      DO ig=1,ngridmx
+       DO l=1,lmax(ig)
+        rpdhd(ig,l)=0.
+        DO k=1,l
+         rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*(pdt_th(ig,k)/zpopsk(ig,k))* &
+     & (zzlev(ig,k+1)-zzlev(ig,k))
+        ENDDO
+        rpdhd(ig,l)=rpdhd(ig,l)+                                        &
+     &    rho(ig,1)*(heatFlux(ig,1)+heatFlux_down(ig,1))
+       ENDDO
+      ENDDO
+      rpdhd(:,nlayermx)=0.
+
+! compute w'teta' from thermals
+
+      wtth(:,:)=rpdhd(:,:)/rho(:,:)
+
+! We get the max heat flux from thermals and add the contribution from the diffusion
+
+      DO ig=1,ngridmx
+        hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:))
+      ENDDO
 ! We follow Spiga et. al 2010 (QJRMS)
 ! ------------
@@ -499,4 +557,11 @@
         call WRITEDIAGFI(ngridmx,'d_u_th',  &
      &         'tendance moment','m/s',1,pdu_th)
+        call WRITEDIAGFI(ngridmx,'wtdif',  &
+     &         'heat flux from diffusion','K.m/s',1,wtdif)
+        call WRITEDIAGFI(ngridmx,'wtth',  &
+     &         'heat flux from thermals','K.m/s',1,wtth)
+        call WRITEDIAGFI(ngridmx,'wttot',  &
+     &         'heat flux PBL','K.m/s',1,wtdif(:,:)+wtth(:,:))
+
       else
 
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 659)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 660)
@@ -353,5 +353,5 @@
       REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
       REAL L_mo(ngridmx),vhf(ngridmx),vvv(ngridmx)
-      REAL zu2(ngridmx)
+      REAL zu2(ngridmx),sensibFlux(ngridmx)
 c=======================================================================
 
@@ -785,5 +785,5 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th
+     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux
 #ifdef MESOSCALE
      &        ,flag_LES
@@ -850,5 +850,5 @@
      $ pplay,pplev,pphi,zpopsk,
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
-     $ dtke_th,hfmax_th,wstar)
+     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
      
          DO l=1,nlayer
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 659)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 660)
@@ -5,5 +5,6 @@
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
-     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax
+     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
+     $                hfmax,sensibFlux
 #ifdef MESOSCALE
      &                ,flag_LES
@@ -143,4 +144,6 @@
       REAL vmr_co2(ngridmx,nlayermx)
       REAL qco2,mmean
+
+      REAL sensibFlux(ngridmx)
 
 #ifdef MESOSCALE
@@ -649,4 +652,7 @@
       pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
 
+      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
+         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig))
+      ENDDO
 
 c-----------------------------------------------------------------------
