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)
      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 "comcstfi.h"
c  des include de la geometrie dynamique pour sorties graphiques
! #include "paramet.h"
! #include "comvert.h"
! #include "comgeom.h"
#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
      REAL, SAVE :: karman=0.4
      INTEGER, SAVE :: icount=0
!$OMP THREADPRIVATE(karman,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)
c     CALL my_25(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,
c    a          pq2,pq2l,zkv,zkh)
        CALL vdif_k(ngrid,nlay,
     s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)

      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