Index: trunk/LMDZ.GENERIC/libf/phygeneric/vdifc_mod.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/vdifc_mod.F	(revision 4036)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/vdifc_mod.F	(revision 4146)
@@ -8,5 +8,5 @@
      &     ptimestep,pcapcal,lecrit,                        
      &     pplay,pplev,pzlay,pzlev,pz0,
-     &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
+     &     pu,pv,pt,ph,pq,ptsrf,pemis,pqsurf,
      &     pdhfi,pdqfi,pfluxsrf,
      &     pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
@@ -18,5 +18,5 @@
       use surfdat_h, only: dryness
       use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
-      use comcstfi_mod, only: g, r, cpp, rcp
+      use comcstfi_mod, only: g, rd_ref, cppd_ref, rcp_ref
       use callkeys_mod, only: water,tracer,nosurf
 
@@ -53,5 +53,6 @@
       REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
       REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
-      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
+      REAL,INTENT(IN) :: pt(ngrid,nlay),ph(ngrid,nlay)
       REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
       REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
@@ -179,9 +180,9 @@
       ENDDO
 
-      zcst1=4.*g*ptimestep/(R*R)
+      zcst1=4.*g*ptimestep/(rd_ref*rd_ref)
       DO ilev=2,nlev-1
          DO ig=1,ngrid
             zb0(ig,ilev)=pplev(ig,ilev)*
-     s           (pplev(ig,1)/pplev(ig,ilev))**rcp /
+     s           (pplev(ig,1)/pplev(ig,ilev))**rcp_ref /
      s           (ph(ig,ilev-1)+ph(ig,ilev))
             zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
@@ -190,5 +191,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
 
@@ -246,6 +247,6 @@
 !     ------------------------------------------------------ 
 
-      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
-     &     ,pu,pv,pq,ph,zcdv_true
+      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,
+     &     pplev,pplay,pu,pv,pt,pq,ph,zcdv_true
      &     ,pq2,zkv,zkh)
 
@@ -405,11 +406,11 @@
          DO ig=1,ngrid
 
-            z1(ig) = pcapcal(ig)*ptsrf(ig) + cpp*zb(ig,1)*zc(ig,1)
-     &           + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
-            z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 
-     &           +zdplanck(ig)
-            ztsrf2(ig) = z1(ig) / z2(ig)
-            pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
-            zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
+           z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zb(ig,1)*zc(ig,1)
+     &          + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
+           z2(ig) = pcapcal(ig) + cppd_ref*zb(ig,1)*(1.-zd(ig,1)) 
+     &          +zdplanck(ig)
+           ztsrf2(ig) = z1(ig) / z2(ig)
+           pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
+           zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
          ENDDO
 
@@ -552,20 +553,20 @@
 ! calculation of h0 and h1
                do ig=1,ngrid
-                  zdq0(ig) = dqsat(ig)
-                  zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
-
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zb(ig,1)*zc(ig,1)
-     &                 + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 
-     &                 +  zb(ig,1)*dryness(ig)*RLVTT*
-     &                 ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
-
-                  z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1))
-     &                 +zdplanck(ig)
-     &                 +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*
-     &                 (1.0-zdq(ig,1))
-
-                  ztsrf2(ig) = z1(ig) / z2(ig)
-                  pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
-                  zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
+                 zdq0(ig) = dqsat(ig)
+                 zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
+
+                 z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zb(ig,1)*zc(ig,1)
+     &                + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 
+     &                +  zb(ig,1)*dryness(ig)*RLVTT*
+     &                ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
+
+                 z2(ig) = pcapcal(ig) + cppd_ref*zb(ig,1)*(1.-zd(ig,1))
+     &                +zdplanck(ig)
+     &                +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*
+     &                (1.0-zdq(ig,1))
+
+                 ztsrf2(ig) = z1(ig) / z2(ig)
+                 pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
+                 zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
                enddo
 
@@ -664,5 +665,5 @@
 
       DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
-	 sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
+	sensibFlux(ig)=cppd_ref*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
       ENDDO