Index: namico_lmdz/simple_physics/phyparam/param/multipl.F
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/param/multipl.F	(revision 4187)
+++ 	(revision )
@@ -1,17 +1,0 @@
-      SUBROUTINE multipl(n,x1,x2,y)
-      IMPLICIT NONE
-c=======================================================================
-c
-c   multiplication de deux vecteurs
-c
-c=======================================================================
-c
-      INTEGER n,i
-      REAL x1(n),x2(n),y(n)
-c
-      DO 10 i=1,n
-	 y(i)=x1(i)*x2(i)
-10    CONTINUE
-c
-      RETURN
-      END
Index: /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F	(revision 4187)
+++ /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F	(revision 4188)
@@ -12,4 +12,5 @@
       USE phys_const
       USE planet
+      USE vdif_mod, ONLY : vdif
 c
       IMPLICIT NONE
@@ -161,5 +162,5 @@
 
 
-      EXTERNAL vdif,convadj
+      EXTERNAL convadj
       EXTERNAL orbite,mucorr
       EXTERNAL ismin,ismax
Index: namico_lmdz/simple_physics/phyparam/param/vdif.F
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/param/vdif.F	(revision 4187)
+++ 	(revision )
@@ -1,299 +1,0 @@
-      subroutiNE vdif(ngrid,nlay,ptime,
-     $                ptimestep,pcapcal,pz0,
-     $                pplay,pplev,pzlay,pzlev,
-     $                pu,pv,ph,ptsrf,pemis,
-     $                pdufi,pdvfi,pdhfi,pfluxsrf,
-     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l,
-     $                lwrite)
-      USE phys_const
-      USE vdif_mod
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Diffusion verticale
-c   Shema implicite
-c   On commence par rajouter au variables x la tendance physique
-c   et on resoult en fait:
-c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
-c
-c   arguments:
-c   ----------
-c
-c   entree:
-c   -------
-c
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-#include "description.h"
-c
-c   arguments:
-c   ----------
-
-      INTEGER ngrid,nlay
-      REAL ptime,ptimestep
-      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
-      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
-      REAL ptsrf(ngrid),pemis(ngrid)
-      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
-      REAL pfluxsrf(ngrid)
-      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
-      REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
-      REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
-      LOGICAL lwrite
-c
-c   local:
-c   ------
-
-      INTEGER ilev,ig,ilay,nlev
-      INTEGER unit,ierr,it1,it2
-      INTEGER cluvdb,putdat,putvdim,setname,setvdim
-      REAL z4st,zdplanck(ngrid),zu2
-      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
-      REAL zcdv(ngrid),zcdh(ngrid)
-      REAL zu(ngrid,nlay),zv(ngrid,nlay)
-      REAL zh(ngrid,nlay)
-      REAL ztsrf2(ngrid)
-      REAL z1(ngrid),z2(ngrid)
-      REAL za(ngrid,nlay),zb(ngrid,nlay)
-      REAL zb0(ngrid,nlay)
-      REAL zc(ngrid,nlay),zd(ngrid,nlay)
-      REAL zcst1
-
-      EXTERNAL coefdifv
-      INTEGER, SAVE :: icount=0
-!$OMP THREADPRIVATE(icount)
-
-c
-c-----------------------------------------------------------------------
-c   initialisations:
-c   ----------------
-
-      nlev=nlay+1
-
-c   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
-c   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
-c   ---------------------------------
-
-      DO ilay=1,nlay
-         DO ig=1,ngrid
-            za(ig,ilay)=
-     s       (pplev(ig,ilay)-pplev(ig,ilay+1))/g
-         ENDDO
-      ENDDO
-
-      zcst1=4.*g*ptimestep/(r*r)
-      DO ilev=2,nlev-1
-         DO ig=1,ngrid
-            zb0(ig,ilev)=pplev(ig,ilev)*
-     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
-     s      (ph(ig,ilev-1)+ph(ig,ilev))
-            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
-     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
-         ENDDO
-      ENDDO
-      DO ig=1,ngrid
-         zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
-      ENDDO
-      IF(lwrite) THEN
-         ig=ngrid/2+1
-         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
-         DO ilay=1,nlay
-            WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
-     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
-         ENDDO
-         PRINT*,'Pression (mbar) ,altitude (km),zb'
-         DO ilev=1,nlay
-            WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
-     s      zb0(ig,ilev)
-         ENDDO
-      ENDIF
-
-c-----------------------------------------------------------------------
-c   2. ajout des tendances physiques:
-c   ------------------------------
-
-      DO ilev=1,nlay
-         DO ig=1,ngrid
-            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
-            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
-            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   3. calcul de  cd :
-c   ----------------
-c
-      CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
-      CALL vdif_k(ngrid,nlay,
-     s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh)
-
-      DO ig=1,ngrid
-         zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
-         zcdv(ig)=zcdv(ig)*sqrt(zu2)
-         zcdh(ig)=zcdh(ig)*sqrt(zu2)
-      ENDDO
-
-      IF(lwrite) THEN
-         PRINT*
-         PRINT*,'Diagnostique diffusion verticale'
-         PRINT*,'coefficients Cd pour v et h'
-         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
-         PRINT*,'coefficients K pour v et h'
-         DO ilev=1,nlay
-            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
-         ENDDO
-      ENDIF
-
-c-----------------------------------------------------------------------
-c   integration verticale pour u:
-c   -----------------------------
-c
-      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
-      CALL multipl(ngrid,zcdv,zb0,zb)
-      DO ig=1,ngrid
-         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
-         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-      ENDDO
-
-      DO ilay=nlay-1,1,-1
-         DO ig=1,ngrid
-            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
-     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-         ENDDO
-      ENDDO
-
-      DO ig=1,ngrid
-         zu(ig,1)=zc(ig,1)
-      ENDDO
-      DO ilay=2,nlay
-         DO ig=1,ngrid
-            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   integration verticale pour v:
-c   -----------------------------
-c
-      DO ig=1,ngrid
-         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
-         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-      ENDDO
-
-      DO ilay=nlay-1,1,-1
-         DO ig=1,ngrid
-            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
-     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-         ENDDO
-      ENDDO
-
-      DO ig=1,ngrid
-         zv(ig,1)=zc(ig,1)
-      ENDDO
-      DO ilay=2,nlay
-         DO ig=1,ngrid
-            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   integration verticale pour h:
-c   -----------------------------
-c
-      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
-      CALL multipl(ngrid,zcdh,zb0,zb)
-      DO ig=1,ngrid
-         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
-         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-      ENDDO
-
-      DO ilay=nlay-1,1,-1
-         DO ig=1,ngrid
-            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
-     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   rajout eventuel de planck dans le shema implicite:
-c   --------------------------------------------------
-
-      z4st=4.*5.67e-8*ptimestep
-c     z4st=0.
-      DO ig=1,ngrid
-         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   calcul le l'evolution de la temperature du sol:
-c   -----------------------------------------------
-
-      DO ig=1,ngrid
-         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
-     s     +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)
-         zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
-         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   integration verticale finale:
-c   -----------------------------
-
-      DO ilay=2,nlay
-         DO ig=1,ngrid
-            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   calcul final des tendances de la diffusion verticale:
-c   -----------------------------------------------------
-
-      DO ilev = 1, nlay
-         DO ig=1,ngrid
-            pdudif(ig,ilev)=(    zu(ig,ilev)-
-     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
-            pdvdif(ig,ilev)=(    zv(ig,ilev)-
-     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
-            pdhdif(ig,ilev)=(    zh(ig,ilev)-
-     $      (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep)    )/ptimestep
-         ENDDO
-      ENDDO
-
-      IF(lwrite) THEN
-         PRINT*
-         PRINT*,'Diagnostique de la diffusion verticale'
-         PRINT*,'h avant et apres diffusion verticale'
-         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
-         DO 3110 ilev=1,nlay
-            PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
-3110     CONTINUE
-      ENDIF
-
-
-c-----------------------------------------------------------------------
-
-      RETURN
-      END
Index: /dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90	(revision 4187)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90	(revision 4188)
@@ -3,11 +3,11 @@
   SAVE
   PRIVATE
-
+  
   REAL, PARAMETER :: karman=0.4
-
-  PUBLIC :: vdif_cd, vdif_k
+  
+  PUBLIC :: vdif
   
 CONTAINS
-
+  
   SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
     !=======================================================================
@@ -96,5 +96,4 @@
     !-----------------------------------------------------------------------
     
-    RETURN
   END SUBROUTINE vdif_cd
   
@@ -149,6 +148,300 @@
     ENDDO
     
-    RETURN
   END SUBROUTINE vdif_k
+
+  
+  SUBROUTINE multipl(n,x1,x2,y)
+!=======================================================================
+!   multiplication de deux vecteurs
+!=======================================================================
+    INTEGER n,i
+    REAL x1(n),x2(n),y(n)
+
+    DO i=1,n
+       y(i)=x1(i)*x2(i)
+    END DO
+  END SUBROUTINE multipl
+
+  SUBROUTINE vdif(ngrid,nlay,ptime, &
+       ptimestep,pcapcal,pz0, &
+       pplay,pplev,pzlay,pzlev, &
+       pu,pv,ph,ptsrf,pemis, &
+       pdufi,pdvfi,pdhfi,pfluxsrf, &
+       pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, &
+       lwrite)
+    USE phys_const
+    
+!=======================================================================
+!
+!   Diffusion verticale
+!   Shema implicite
+!   On commence par rajouter au variables x la tendance physique
+!   et on resoult en fait:
+!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
+!
+!   arguments:
+!   ----------
+!
+!   entree:
+!   -------
+!
+!
+!=======================================================================
+
+!
+!   arguments:
+!   ----------
+
+    INTEGER ngrid,nlay
+    REAL ptime,ptimestep
+    REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+    REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+    REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+    REAL ptsrf(ngrid),pemis(ngrid)
+    REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
+    REAL pfluxsrf(ngrid)
+    REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
+    REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
+    REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
+    LOGICAL lwrite
+    !
+    !   local:
+    !   ------
+    
+    INTEGER ilev,ig,ilay,nlev
+    INTEGER unit,ierr,it1,it2
+    INTEGER cluvdb,putdat,putvdim,setname,setvdim
+    REAL z4st,zdplanck(ngrid),zu2
+    REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
+    REAL zcdv(ngrid),zcdh(ngrid)
+    REAL zu(ngrid,nlay),zv(ngrid,nlay)
+    REAL zh(ngrid,nlay)
+    REAL ztsrf2(ngrid)
+    REAL z1(ngrid),z2(ngrid)
+    REAL za(ngrid,nlay),zb(ngrid,nlay)
+    REAL zb0(ngrid,nlay)
+    REAL zc(ngrid,nlay),zd(ngrid,nlay)
+    REAL zcst1
+    
+    !
+    !-----------------------------------------------------------------------
+    !   initialisations:
+    !   ----------------
+    
+    nlev=nlay+1
+    
+    !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
+    !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
+    !   ---------------------------------
+    
+    DO ilay=1,nlay
+       DO ig=1,ngrid
+          za(ig,ilay) = (pplev(ig,ilay)-pplev(ig,ilay+1))/g
+       ENDDO
+    ENDDO
+    
+    zcst1=4.*g*ptimestep/(r*r)
+    DO ilev=2,nlev-1
+       DO ig=1,ngrid
+          zb0(ig,ilev)=pplev(ig,ilev) &
+               *(pplev(ig,1)/pplev(ig,ilev))**rcp &
+               /(ph(ig,ilev-1)+ph(ig,ilev))
+          zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev) &
+               / (pplay(ig,ilev-1)-pplay(ig,ilev))
+       ENDDO
+    ENDDO
+    DO ig=1,ngrid
+       zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
+    ENDDO
+    IF(lwrite) THEN
+       ig=ngrid/2+1
+       PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
+       DO ilay=1,nlay
+          WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), &
+               pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
+       ENDDO
+       PRINT*,'Pression (mbar) ,altitude (km),zb'
+       DO ilev=1,nlay
+          WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), &
+               zb0(ig,ilev)
+       ENDDO
+    ENDIF
+    
+    !-----------------------------------------------------------------------
+    !   2. ajout des tendances physiques:
+    !   ------------------------------
+    
+    DO ilev=1,nlay
+       DO ig=1,ngrid
+          zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
+          zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+          zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   3. calcul de  cd :
+    !   ----------------
+    !
+    CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
+    CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh)
+    
+    DO ig=1,ngrid
+       zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
+       zcdv(ig)=zcdv(ig)*sqrt(zu2)
+       zcdh(ig)=zcdh(ig)*sqrt(zu2)
+    ENDDO
+    
+    IF(lwrite) THEN
+       PRINT*
+       PRINT*,'Diagnostique diffusion verticale'
+       PRINT*,'coefficients Cd pour v et h'
+       PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
+       PRINT*,'coefficients K pour v et h'
+       DO ilev=1,nlay
+          PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
+       ENDDO
+    ENDIF
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour u:
+    !   -----------------------------
+    
+    CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
+    CALL multipl(ngrid,zcdv,zb0,zb)
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig) = 1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay) = (za(ig,ilay)*zu(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay) = zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    DO ig=1,ngrid
+       zu(ig,1)=zc(ig,1)
+    ENDDO
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour v:
+    !   -----------------------------
+    !
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    DO ig=1,ngrid
+       zv(ig,1)=zc(ig,1)
+    ENDDO
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour h:
+    !   -----------------------------
+    !
+    CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
+    CALL multipl(ngrid,zcdh,zb0,zb)
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   rajout eventuel de planck dans le shema implicite:
+    !   --------------------------------------------------
+    
+    z4st=4.*5.67e-8*ptimestep
+    DO ig=1,ngrid
+       zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   calcul le l'evolution de la temperature du sol:
+    !   -----------------------------------------------
+    
+    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)
+       zh(ig,1) = zc(ig,1)+zd(ig,1)*ztsrf2(ig)
+       pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale finale:
+    !   -----------------------------
+    
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   calcul final des tendances de la diffusion verticale:
+    !   -----------------------------------------------------
+    
+    DO ilev = 1, nlay
+       DO ig=1,ngrid
+          pdudif(ig,ilev) = ( zu(ig,ilev) &
+               - (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
+          pdvdif(ig,ilev) = ( zv(ig,ilev) &
+               - (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
+          pdhdif(ig,ilev) = ( zh(ig,ilev) &
+               - (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep
+       ENDDO
+    ENDDO
+    
+    IF(lwrite) THEN
+       PRINT*
+       PRINT*,'Diagnostique de la diffusion verticale'
+       PRINT*,'h avant et apres diffusion verticale'
+       PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
+       DO  ilev=1,nlay
+          PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
+       ENDDO
+    END IF
+
+  END SUBROUTINE vdif
   
 END MODULE vdif_mod
