subroutine thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,masse  &
     &    ,fraca,larga  &
     &    ,u,v,du,dv,ua,va,lev_out)
      implicit none

!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
!=======================================================================


      integer ngrid,nlay

      real ptimestep
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real fraca(ngrid,nlay+1)
      real larga(ngrid)
      real entr(ngrid,nlay)
      real u(ngrid,nlay)
      real ua(ngrid,nlay)
      real du(ngrid,nlay)
      real v(ngrid,nlay)
      real va(ngrid,nlay)
      real dv(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2
      real wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
      real gamma0,gamma(ngrid,nlay+1)
      real ue(ngrid,nlay),ve(ngrid,nlay)
      real dua,dva
      integer iter

      integer ig,k

!   calcul du detrainement

      do k=1,nlay
         do ig=1,ngrid
            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
         enddo
      enddo

!   calcul de la valeur dans les ascendances
      do ig=1,ngrid
         ua(ig,1)=u(ig,1)
         va(ig,1)=v(ig,1)
         ue(ig,1)=u(ig,1)
         ve(ig,1)=v(ig,1)
      enddo

      do k=2,nlay
         do ig=1,ngrid
            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     &         1.e-5*masse(ig,k)) then
!   On itère sur la valeur du coeff de freinage.
!              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
               gamma0=masse(ig,k)  &
     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
     &         *0.5/larga(ig)  &
     &         *1.
!    s         *0.5
!              gamma0=0.
               zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
               zf=0.
               zf2=1./(1.-zf)
!   la première fois on multiplie le coefficient de freinage
!   par le module du vent dans la couche en dessous.
               dua=ua(ig,k-1)-u(ig,k-1)
               dva=va(ig,k-1)-v(ig,k-1)
               do iter=1,5
!   On choisit une relaxation lineaire.
                  gamma(ig,k)=gamma0
!   On choisit une relaxation quadratique.
                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  &
     &               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))  &
     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
     &                 +gamma(ig,k))
                  va(ig,k)=(fm(ig,k)*va(ig,k-1)  &
     &               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))  &
     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
     &                 +gamma(ig,k))
!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
                  dua=ua(ig,k)-u(ig,k)
                  dva=va(ig,k)-v(ig,k)
                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
               enddo
            else
               ua(ig,k)=u(ig,k)
               va(ig,k)=v(ig,k)
               ue(ig,k)=u(ig,k)
               ve(ig,k)=v(ig,k)
               gamma(ig,k)=0.
            endif
         enddo
      enddo

      do k=2,nlay
         do ig=1,ngrid
            wud(ig,k)=fm(ig,k)*ue(ig,k)
            wvd(ig,k)=fm(ig,k)*ve(ig,k)
         enddo
      enddo
      do ig=1,ngrid
         wud(ig,1)=0.
         wud(ig,nlay+1)=0.
         wvd(ig,1)=0.
         wvd(ig,nlay+1)=0.
      enddo

      do k=1,nlay
         do ig=1,ngrid
            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
     &               -wud(ig,k)+wud(ig,k+1))  &
     &               /masse(ig,k)
            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)  &
     &               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)  &
     &               -wvd(ig,k)+wvd(ig,k+1))  &
     &               /masse(ig,k)
         enddo
      enddo

      return
      end