--- trunk/libf/phylmd/thermcell.f	2008/02/27 13:16:39	3
+++ trunk/libf/phylmd/thermcell.f	2010/04/06 17:52:58	32
@@ -45,7 +45,7 @@
       REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
       REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
       REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
-      REAL pplay(ngrid,nlay)
+      REAL, intent(in):: pplay(ngrid,nlay)
       real, intent(in):: pplev(ngrid,nlay+1)
       real pphi(ngrid,nlay)
 
@@ -722,14 +722,6 @@
             endif
             if(.not.masse(ig,l).ge.1.e-10
      s         .or..not.masse(ig,l).le.1.e4) then
-c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
-c    s         ,'   M=',masse(ig,l)
-c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
-c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
-c     print*,'zlev(ig,l+1),zlev(ig,l)'
-c    s                ,zlev(ig,l+1),zlev(ig,l)
-c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
-c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
             endif
             if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
 c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
@@ -933,120 +925,7 @@
 
       return
       end
-      subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
-     .    ,fraca,larga
-     .    ,u,v,du,dv,ua,va)
-      use dimens_m
-      use dimphy
-      implicit none
-
-c=======================================================================
-c
-c   Calcul du transport verticale dans la couche limite en presence
-c   de "thermiques" explicitement representes
-c   calcul du dq/dt une fois qu'on connait les ascendances
-c
-c=======================================================================
-
-
-      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)
-
-      real qa(klon,klev),detr(klon,klev)
-      real wvd(klon,klev+1),wud(klon,klev+1)
-      real gamma0,gamma(klon,klev+1)
-      real dua,dva
-      integer iter
-
-      integer ig,k
-
-c   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
-
-c   calcul de la valeur dans les ascendances
-      do ig=1,ngrid
-         ua(ig,1)=u(ig,1)
-         va(ig,1)=v(ig,1)
-      enddo
-
-      do k=2,nlay
-         do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
-     s         1.e-5*masse(ig,k)) then
-c   On itère sur la valeur du coeff de freinage.
-c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
-               gamma0=masse(ig,k)
-     s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
-     s         *0.5/larga(ig)
-c              gamma0=0.
-c   la première fois on multiplie le coefficient de freinage
-c   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
-                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
-                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
-     s               +(entr(ig,k)+gamma(ig,k))*u(ig,k))
-     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
-                  va(ig,k)=(fm(ig,k)*va(ig,k-1)
-     s               +(entr(ig,k)+gamma(ig,k))*v(ig,k))
-     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
-c                 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)
-               enddo
-            else
-               ua(ig,k)=u(ig,k)
-               va(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)*u(ig,k)
-            wvd(ig,k)=fm(ig,k)*v(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)
-     s               -(entr(ig,k)+gamma(ig,k))*u(ig,k)
-     s               -wud(ig,k)+wud(ig,k+1))
-     s               /masse(ig,k)
-            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
-     s               -(entr(ig,k)+gamma(ig,k))*v(ig,k)
-     s               -wvd(ig,k)+wvd(ig,k+1))
-     s               /masse(ig,k)
-         enddo
-      enddo
-
-      return
-      end
       subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
      .    ,q,dq,qa)
       use dimens_m