Index: trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90	(revision 4036)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90	(revision 4146)
@@ -8,5 +8,5 @@
           ptimestep,pcapcal,                    &   
           pplay,pplev,pzlay,pzlev,pz0,                 &
-          pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
+          pu,pv,pt,ph,ppopsk,pq,ptsrf,pemis,pqsurf,       &
           pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
@@ -17,5 +17,5 @@
       use surfdat_h, only: dryness
       use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
-      use comcstfi_mod, only: rcp, g, r, cpp
+      use comcstfi_mod, only: g, rd_ref, cppd_ref
       use callkeys_mod, only: water,tracer,nosurf,kmixmin
       use turb_mod, only : ustar
@@ -61,5 +61,5 @@
       REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
       REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
-      REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
+      REAL,INTENT(IN) :: pt(ngrid,nlay),ph(ngrid,nlay),ppopsk(ngrid,nlay)
       REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
       REAL,INTENT(IN) :: pemis(ngrid)
@@ -173,10 +173,10 @@
          DO ig=1,ngrid
             zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
-	    zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
+	    zExner(ig,ilay)=ppopsk(ig,ilay)
 	    zovExner(ig,ilay)=1./ppopsk(ig,ilay)
          ENDDO
       ENDDO
 
-      zcst1=4.*g*ptimestep/(R*R)
+      zcst1=4.*g*ptimestep/(rd_ref*rd_ref)
       DO ilev=2,nlev-1
          DO ig=1,ngrid
@@ -186,5 +186,5 @@
       ENDDO
       DO ig=1,ngrid
-         zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig))
+         zb0(ig,1)=ptimestep*pplev(ig,1)/(rd_ref*ptsrf(ig))
       ENDDO
       dqsdif_total(:)=0.0
@@ -199,5 +199,5 @@
             zv(ig,ilev)=pv(ig,ilev)
             zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
-            zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
+            zh(ig,ilev)=ph(ig,ilev)!+pdtfi(ig,ilev)*ptimestep*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
         ENDDO
       ENDDO
@@ -248,6 +248,7 @@
 !     ------------------------------------------------------ 
 
-      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pu,pv,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
-      
+      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pplev,pplay, & 
+                   pu,pv,zt,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated
+                       
 !     Adding eddy mixing to mimic 3D general circulation in 1D
 !     R. Wordsworth & F. Forget (2010)
@@ -273,5 +274,4 @@
 !               zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev)) 
 !            end do
-
             do ig=1,ngrid
                zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
@@ -436,7 +436,7 @@
 
          DO ig=1,ngrid
-            z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
+            z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
                 + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig) 
-            z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) 
+            z2(ig) = pcapcal(ig)+zdplanck(ig)+cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) 
             ztsrf(ig) = z1(ig) / z2(ig)
             pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
@@ -573,9 +573,9 @@
                   zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
 
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
+                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
 		      + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
                       + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
 
-                  z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
+                  z2(ig) = pcapcal(ig) + cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
                       + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1))
 
@@ -598,8 +598,8 @@
 
                       !recompute surface temperature  
-                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
+                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
 		        + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
                         + RLVTT*dqsdif_total(ig)
-                      z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
+                      z2(ig) = pcapcal(ig) + cppd_ref*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
                         + zdplanck(ig)
                       ztsrf(ig) = z1(ig) / z2(ig)
@@ -720,5 +720,5 @@
       
       DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
-	 sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
+	 sensibFlux(ig)=cppd_ref*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
       ENDDO