Index: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 648)
+++ trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 652)
@@ -3,5 +3,5 @@
 !
       SUBROUTINE calltherm_interface (firstcall, &
-     & long,lati,zzlev,zzlay, &
+     & zzlev,zzlay, &
      & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
      & pplay,pplev,pphi,zpopsk, &
@@ -21,4 +21,5 @@
 !--------------------------------------------------------
 
+!      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
       REAL, INTENT(IN) :: ptimestep
       REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
@@ -32,5 +33,4 @@
       REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
       REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
-      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
       REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
 
@@ -363,5 +363,5 @@
       if (dtke_thermals) then
       detrmod(:,:)=0.
-      ndt=1
+      ndt=10
       do l=1,zlmax
          do ig=1,ngridmx
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 648)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 652)
@@ -846,5 +846,5 @@
  
         call calltherm_interface(firstcall,
-     $ long,lati,zzlev,zzlay,
+     $ zzlev,zzlay,
      $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
      $ pplay,pplev,pphi,zpopsk,
Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 648)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 652)
@@ -857,5 +857,6 @@
             if (l.le.lmax(ig)) then
                 if (zw2(ig,l).lt.0.)then
-                  print*,'pb2 zw2<0'
+!                  print*,'pb2 zw2<0',zw2(ig,l)
+                  zw2(ig,l)=0.
                 endif
                 zw2(ig,l)=sqrt(zw2(ig,l))
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 648)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 652)
@@ -384,5 +384,5 @@
       CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
      s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ust
-     s   ,8)
+     s   ,9)
 
       ENDIF
Index: trunk/LMDZ.MARS/libf/phymars/yamada4.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 648)
+++ trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 652)
@@ -120,4 +120,8 @@
       REAL teta(ngrid,nlay)
       REAL pq(ngrid,nlay,nqmx)
+      REAL kminfact
+      INTEGER i
+      REAL ztimestep
+      INTEGER, SAVE :: ndt
 
       nlev=nlay+1
@@ -146,4 +150,5 @@
         endif
       allocate(l0(ngrid))
+      ndt=ceiling(3840./(3699.*24./dt))
       firstcall=.false.
       endif !of if firstcall
@@ -153,14 +158,14 @@
 c.......................................................................
 
-      if (ico2.ne.0) then
-!     Special case if one of the tracers is CO2 gas
-         DO k=1,nlay
-           DO ig=1,ngrid
-            teta(ig,k) = phc(ig,k)*(A*pq(ig,k,ico2)+B)
-           ENDDO
-         ENDDO
-       else
+!      if (ico2.ne.0) then
+!!     Special case if one of the tracers is CO2 gas
+!         DO k=1,nlay
+!           DO ig=1,ngrid
+!            teta(ig,k) = phc(ig,k)*(A*pq(ig,k,ico2)+B)
+!           ENDDO
+!         ENDDO
+!       else
           teta(:,:)=phc(:,:)
-       end if
+!       end if
       
       if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
@@ -402,8 +407,10 @@
 !====================================================================
 
+      ztimestep=dt/real(ndt)
+      do i=1,ndt
 
 !  Calcul de l,  km, au pas precedent
       do k=2,nlay
-                                                          do ig=1,ngrid
+       do ig=1,ngrid
 !        print*,'SMML=',sm(ig,k),l(ig,k)
          delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
@@ -418,5 +425,5 @@
      s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
 ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
-         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
+         aa(ig,k)=aa1*ztimestep/(delta(ig,k)*l(ig,k))
 !     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
          qpre=sqrt(q2(ig,k))
@@ -440,5 +447,5 @@
 
 !     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
-                                                          enddo
+       enddo
       enddo
 
@@ -446,4 +453,24 @@
       q2(:,nlay+1)=q2(:,nlay)
 
+      if (iflag_pbl .eq. 9) then
+      do k=2,nlay
+      do ig=1,ngrid
+        zq=sqrt(q2(ig,k))
+        km(ig,k)=l(ig,k)*zq*sm(ig,k)
+        kn(ig,k)=km(ig,k)*alpha(ig,k)
+        kq(ig,k)=l(ig,k)*zq*0.2
+      enddo
+      enddo
+      km(:,nlay+1)=km(:,nlay)
+      kn(:,nlay+1)=kn(:,nlay)
+      kq(:,nlay+1)=kq(:,nlay)
+
+      q2(:,1)=q2(:,2)
+       call vdif_q2(ztimestep,g,rconst,ngrid,nlay,plev,temp,kq,q2)
+
+      endif ! of if iflag_pbl eq 9
+
+      enddo !of i=1,ndt
+
       endif ! Fin du cas 8
 
@@ -453,4 +480,5 @@
 !   Calcul des coefficients de melange
 !====================================================================
+      if (iflag_pbl .ne. 9) then
       do k=2,nlay
 !     print*,'k=',k
@@ -471,8 +499,10 @@
 
 ! Transport diffusif vertical de la TKE.
-      if (iflag_pbl.ge.9) then
-!       print*,'YAMADA VDIF'
-        q2(:,1)=q2(:,2)
-        call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
+!      if (iflag_pbl.ge.9) then
+!!       print*,'YAMADA VDIF'
+!        q2(:,1)=q2(:,2)
+!        call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
+!      endif
+
       endif
 
@@ -484,24 +514,34 @@
 !   D'apres Holtslag Boville.
 
-! MARS : deactivating pblhmin following F.H. concerns
-
-!                                                          do ig=1,ngrid
+! MARS
+!       callkmin=.true.
+!       call getin("callkmin",callkmin)
+!       IF (callkmin) THEN
+                                                          do ig=1,ngrid
 !      coriol(ig)=1.e-4
 !      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
-!                                                          enddo
-!
+
+       if (ngrid .eq. 1) then
+       kminfact=0.3
+       else
+       kminfact=0.45
+       endif
+
+       pblhmin(ig)=kminfact*0.07*MAX(ustar(ig),1.e-3)/1.e-4
+                                                   enddo
 !      print*,'pblhmin ',pblhmin
 !CTest a remettre 21 11 02
 ! test abd 13 05 02      if(0.eq.1) then
 !      if(0.eq.1) then
-!      do k=2,nlay
-!         do ig=1,ngrid
-!            if (teta(ig,2).gt.teta(ig,1)) then
-!               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
+      do k=2,nlay
+         do ig=1,ngrid
+            if (teta(ig,2).gt.teta(ig,1)) then
+               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
 !               kmin=kap*zlev(ig,k)*qmin
-!            else
-!               kmin=-1. ! kmin n'est utilise que pour les SL stables.
-!            endif 
-!            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
+               kmin=fl(zlev(ig,k),l0(ig),qmin**2,n2(ig,k))*qmin
+            else
+               kmin=-1. ! kmin n'est utilise que pour les SL stables.
+            endif
+            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
 !               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
 !     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
@@ -509,11 +549,19 @@
 !               km(ig,k)=kmin
 !               kq(ig,k)=kmin
+
+               kn(ig,k)=kmin*alpha(ig,k)
+               km(ig,k)=kmin
+               kq(ig,k)=kmin*0.2
 !   la longueur de melange est suposee etre l= kap z
 !   K=l q Sm d'ou q2=(K/l Sm)**2
 !               q2(ig,k)=(qmin/sm(ig,k))**2
-!            endif
-!         enddo
-!      enddo
+               q2(ig,k)=(kmin/
+     &     (fl(zlev(ig,k),l0(ig),qmin**2,n2(ig,k))*sm(ig,k)))**2
+            endif
+         enddo
+      enddo
 !      endif
+
+!      ENDIF
 
 !   Diagnostique pour stokage
