Index: /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 528)
+++ /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 529)
@@ -1,5 +1,5 @@
-      SUBROUTINE pbl_parameters(ngrid,nlay,pz0, 
-     & pg,pz,pu,pv,wmax,hfmax,zmax,pts,ph,z_out,
-     & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
+      SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, 
+     & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out,
+     & Teta_out,u_out,ustar,tstar,L_mo,vhf,vvv)
       IMPLICIT NONE
 !=======================================================================
@@ -21,9 +21,9 @@
 !     pz0(ngrid)       surface roughness length
 !     pg               gravity (m s -2)
-!     pz(ngrid,nlay)   height of layers
+!     zzlay(ngrid,nlay)   height of mid-layers
+!     zzlev(ngrid,nlay+1)   height of layers interfaces
 !     pu(ngrid,nlay)   u component of the wind
 !     pv(ngrid,nlay)   v component of the wind
-!     wmax(ngrid)      maximum vertical velocity in thermals (might not be needed
-!                                          if the computation of w* works)
+!     wstar_in(ngrid)  free convection velocity in thermals
 !     hfmax(ngrid)     maximum vertical turbulent heat flux in thermals
 !     zmax(ngrid)      height reached by the thermals (pbl height)
@@ -56,12 +56,18 @@
 
       INTEGER, INTENT(IN) :: ngrid,nlay
-      REAL, INTENT(IN) :: pz0(ngrid)
-      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
+      REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay)
+      REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay)
       REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
-      REAL, INTENT(IN) :: wmax(ngrid),hfmax(ngrid),zmax(ngrid)
+      REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid)
       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
       REAL, INTENT(IN) :: z_out
+
+!    Outputs:
+!    --------
+
       REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)
-      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid),wstar(ngrid)
+      REAL T_out(ngrid)
+      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid)
+      REAL wstar(ngrid)
       REAL, INTENT(OUT) :: L_mo(ngrid)
 
@@ -139,5 +145,6 @@
        ite=0.
        residual=abs(pz0tcomp(ig)-pz0t)
-       zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1),wmax(ig)**2)
+       zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
+     &                                ,(0.3*wstar_in(ig))**2)
 
        DO WHILE((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))
@@ -147,6 +154,6 @@
             ! Richardson number formulation proposed by D.E. England et al. (1995)
             rib(ig) = (pg/ph(ig,1))
-     &        *sqrt(pz(ig,1)*pz0(ig))
-     &        *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
+     &        *sqrt(zzlev(ig,2)*pz0(ig))
+     &        *(((log(zzlev(ig,2)/pz0(ig)))**2)/(log(zzlev(ig,2)/pz0t)))
      &        *(ph(ig,1)-pts(ig))/zu2(ig)
          ELSE
@@ -156,6 +163,6 @@
          ENDIF
 
-         z1z0=pz(ig,1)/pz0(ig)
-         z1z0t=pz(ig,1)/pz0t
+         z1z0=zzlev(ig,2)/pz0(ig)
+         z1z0t=zzlev(ig,2)/pz0t
 
          cdn(ig)=karman/log(z1z0)
@@ -310,4 +317,10 @@
       ENDIF
 
+              T_out(:) = Teta_out(:)*(exp(
+     &   (zout/zzlay(:,1))*(log(pplay(:,1)/ps))
+     &                             )
+     &                         )**rcp
+
+
 !------------------------------------------------------------------------
 !------------------------------------------------------------------------
@@ -325,6 +338,6 @@
       DO k=1,nlay-1
          DO ig=1, ngrid
-            IF (abs(zmax(ig)-pz(ig,k)) .lt. 
-     &                  abs(zmax(ig)-pz(ig,pbl_height_index(ig)))) THEN
+            IF (abs(zmax(ig)-zzlay(ig,k)) .lt. 
+     &              abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN
                pbl_height_index(ig)=k
             ENDIF
@@ -337,5 +350,5 @@
       DO k=1,nlay-1
          DO ig=1, ngrid
-            dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(pz(ig,k+1)-pz(ig,k))
+         dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k))
          ENDDO
       ENDDO
@@ -424,13 +437,13 @@
      &              'Monin Obukhov length','m',
      &                         dimout,L_mo)
-         call WRITEDIAGFI(ngrid,'w_star',
-     &              'Free convection velocity','m',
-     &                         dimout,wstar)
-         call WRITEDIAGFI(ngrid,'z_i',
-     &              'PBL height','m',
-     &                         dimout,zmax)
-         call WRITEDIAGFI(ngrid,'hf_max',
-     &              'Maximum vertical heat flux','m',
-     &                         dimout,hfmax)
+!         call WRITEDIAGFI(ngrid,'w_star',
+!     &              'Free convection velocity','m',
+!     &                         dimout,wstar)
+!         call WRITEDIAGFI(ngrid,'z_i',
+!     &              'PBL height','m',
+!     &                         dimout,zmax)
+!         call WRITEDIAGFI(ngrid,'hf_max',
+!     &              'Maximum vertical heat flux','m',
+!     &                         dimout,hfmax)
          call WRITEDIAGFI(ngrid,'vvv',
      &              'Vertical velocity variance at zout','m',
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 528)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 529)
@@ -777,5 +777,10 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th)
+     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th
+#ifdef MESOSCALE
+     &        ,flag_LES
+#endif
+     &        )
+
 
 #ifdef MESOSCALE
@@ -1848,6 +1853,6 @@
          if (calltherm .and. (z_out .gt. 0.)) then
 
-         call pbl_parameters(ngrid,nlayer,z0,
-     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+         call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
+     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
      & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
 
@@ -1938,9 +1943,9 @@
 ! THERMALS STUFF 1D
 
-         z_out=0.
+         z_out=1.
          if (calltherm .and. (z_out .gt. 0.)) then
 
-         call pbl_parameters(ngrid,nlayer,z0,
-     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+         call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
+     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
      & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
 
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 528)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 529)
@@ -5,5 +5,10 @@
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
-     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax)
+     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax
+#ifdef MESOSCALE
+     &                ,flag_LES
+#endif
+     &                )
+       USE ioipsl_getincom
       IMPLICIT NONE
 
@@ -137,4 +142,7 @@
       REAL qco2,mmean
 
+#ifdef MESOSCALE
+      LOGICAL flag_LES     !! pour LES avec isfflx!=0
+#endif
 c    ** un petit test de coherence
 c       --------------------------
@@ -242,5 +250,6 @@
          DO ilev=1,nlay
             DO ig=1,ngrid
-              qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep
+              qco2=MAX(1.E-30
+     &             ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
 c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
               mmean=1/(A*qco2 +B)
@@ -339,8 +348,7 @@
         ENDIF
 
-
 ! Some usefull diagnostics for the new surface layer parametrization :
- 
-!         call WRITEDIAGFI(ngridmx,'pcdv',
+
+!        call WRITEDIAGFI(ngridmx,'pcdv',
 !     &              'momentum drag','no units',
 !     &                         2,zcdv_true)
@@ -361,9 +369,9 @@
 c    ** schema de diffusion turbulente dans la couche limite
 c       ---------------------------------------------------- 
-!
+
        CALL vdif_kc(ptimestep,g,pzlev,pzlay
      &              ,pu,pv,ph,zcdv_true
      &              ,pq2,zkv,zkh,zq)
-!
+
 !      pt(:,:)=ph(:,:)*ppopsk(:,:)
 !      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
@@ -519,5 +527,4 @@
       ENDDO
 
-
 ! calcul de zc et zd pour la couche top en prenant en compte le terme
 ! de variation de masse (on fait une boucle pour que ça converge)
@@ -526,5 +533,7 @@
 
       llnt(:)=1
-
+#ifdef MESOSCALE
+      IF (.not.flag_LES) THEN
+#endif
       DO ig=1,ngrid
          DO l=1,nlay
@@ -537,4 +546,8 @@
 5     continue
       ENDDO
+
+#ifdef MESOSCALE
+      ENDIF
+#endif
 
       DO ig=1,ngrid
@@ -589,5 +602,4 @@
          ztsrf2(ig)=z1(ig)/z2(ig)
 !         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep  !incremented outside loop
-
             zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
 
@@ -633,9 +645,14 @@
 c     Using the wind modified by friction for lifting and  sublimation
 c     ----------------------------------------------------------------
-c        DO ig=1,ngrid
-c          zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
-c          zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
-c          zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
-c        ENDDO
+
+!     This is computed above and takes into account surface-atmosphere flux 
+!     enhancement by subgrid gustiness and atmospheric-stability related
+!     variations of transfer coefficients.
+
+!        DO ig=1,ngrid
+!          zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
+!          zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
+!          zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
+!        ENDDO
 
 c       Calcul du flux vertical au bas de la premiere couche (dust) :
