Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 498)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 499)
@@ -6,5 +6,5 @@
      & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
      & pplay,pplev,pphi,zpopsk, &
-     & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wmax)
+     & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar)
 
        USE ioipsl_getincom
@@ -45,4 +45,5 @@
       REAL, INTENT(OUT) :: zmaxth(ngridmx)
       REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1)
+      REAL, INTENT(OUT) :: wstar(ngridmx)
 
 !--------------------------------------------------------
@@ -66,5 +67,5 @@
       REAL zdz(ngridmx,nlayermx)
       LOGICAL qtransport_thermals,dtke_thermals
-      INTEGER l,ig,iq
+      INTEGER l,ig,iq,ii(1)
       CHARACTER (LEN=20) :: modname
 
@@ -97,4 +98,5 @@
       REAL buoyancyEst(ngridmx,nlayermx)
       REAL hfmax(ngridmx),wmax(ngridmx)
+      REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
 
 !---------------------------------------------------------
@@ -352,4 +354,40 @@
            pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
 
+
+! **********************************************************************
+! Compute the free convection velocity scale for vdifc
+! **********************************************************************
+
+
+! Potential temperature gradient
+
+      dteta(:,nlayermx)=0.
+      DO l=1,nlayermx-1
+         DO ig=1, ngridmx
+            dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l))          &
+     &              /(zzlay(ig,l+1)-zzlay(ig,l))
+         ENDDO
+      ENDDO
+
+! Computation of the pbl mixed layer temperature
+
+      DO ig=1, ngridmx
+         ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
+         pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
+      ENDDO
+
+! We follow Spiga et. al 2010 (QJRMS)
+! ------------
+
+      DO ig=1, ngridmx
+         IF (zmax(ig) .gt. 0.) THEN
+            wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.)
+         ELSE
+            wstar(ig)=0.
+         ENDIF
+      ENDDO
+
+
+
 ! **********************************************************************
 ! Diagnostics
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 498)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 499)
@@ -332,5 +332,5 @@
 
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
-      REAL, SAVE :: wmax_th(ngridmx)
+      REAL, SAVE :: wstar(ngridmx)
       REAL, SAVE :: hfmax_th(ngridmx)
       REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
@@ -342,5 +342,5 @@
       REAL z_out                          ! height of interpolation between z0 and z1 [meters]
       REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
-      REAL L_mo(ngridmx),wstar(ngridmx),vhf(ngridmx),vvv(ngridmx)
+      REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx)
       REAL zu2(ngridmx)
 c=======================================================================
@@ -359,5 +359,5 @@
          fluxrad(:)=0
 
-         wmax_th(:)=0.
+         wstar(:)=0.
 
 c        read startfi 
@@ -747,8 +747,8 @@
           DO ig=1, ngridmx
              IF (zh(ig,1) .lt. tsurf(ig)) THEN
-               wmax_th(ig)=1.
+               wstar(ig)=1.
                hfmax_th(ig)=0.2
              ELSE
-               wmax_th(ig)=0.
+               wstar(ig)=0.
                hfmax_th(ig)=0.
              ENDIF      
@@ -768,5 +768,5 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh,hfmax_th)
+     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th)
 
 #ifdef MESOSCALE
@@ -829,5 +829,5 @@
      $ pplay,pplev,pphi,zpopsk,
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
-     $ dtke_th,hfmax_th,wmax_th)
+     $ dtke_th,hfmax_th,wstar)
  
          DO l=1,nlayer
@@ -858,5 +858,5 @@
         else   !of if calltherm
         lmax_th(:)=0
-        wmax_th(:)=0.
+        wstar(:)=0.
         hfmax_th(:)=0.
         lmax_th_out(:)=0.
@@ -1812,6 +1812,6 @@
 
          call pbl_parameters(ngrid,nlayer,z0,
-     & g,zzlay,zu,zv,wmax_th,hfmax_th,zmax_th,tsurf,zh,z_out,
-     & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
+     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+     & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
 
          else
@@ -1851,7 +1851,7 @@
      &              'maximum TH heat flux','K.m/s',
      &                         2,hfmax_th)
-        call WRITEDIAGFI(ngridmx,'wmax_th',
+        call WRITEDIAGFI(ngridmx,'wstar',
      &              'maximum TH vertical velocity','m/s',
-     &                         2,wmax_th)
+     &                         2,wstar)
 
          endif
@@ -1905,6 +1905,6 @@
 
          call pbl_parameters(ngrid,nlayer,z0,
-     & g,zzlay,zu,zv,wmax_th,hfmax_th,zmax_th,tsurf,zh,z_out,
-     & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
+     & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
+     & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
 
          else
@@ -1927,7 +1927,7 @@
      &              'maximum TH heat flux','K.m/s',
      &                         0,hfmax_th)
-        call WRITEDIAGFI(ngridmx,'wmax_th',
+        call WRITEDIAGFI(ngridmx,'wstar',
      &              'maximum TH vertical velocity','m/s',
-     &                         0,wmax_th)
+     &                         0,wstar)
 
          co2col(:)=0.
Index: /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 498)
+++ /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 499)
@@ -1,4 +1,4 @@
       SUBROUTINE vdif_cd(ngrid,nlay,pz0,
-     &           pg,pz,pu,pv,wmax,pts,ph,pcdv,pcdh)
+     &           pg,pz,pu,pv,wstar,pts,ph,pcdv,pcdh)
       IMPLICIT NONE
 c=======================================================================
@@ -46,5 +46,5 @@
       REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
-      REAL, INTENT(IN) :: wmax(ngrid)
+      REAL, INTENT(IN) :: wstar(ngrid)
       REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient
 
@@ -146,10 +146,11 @@
 c Richardson number formulation proposed by D.E. England et al. (1995)
 
-!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wmax(ig)**2)
+!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wstar(ig)**2)
 !         zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
-         zu2(ig)=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)
-!       zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.5*wmax(ig))**2
-
-               ! we add the wmax to simulate
+         zu2(ig)=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),             &
+     &      (0.3*wstar(ig))**2)
+!       zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.5*wstar(ig))**2
+
+               ! we add the wstar to simulate
                ! bulk Ri changes due to subgrid wind feeding the thermals
 
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 498)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 499)
@@ -5,5 +5,5 @@
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
-     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true,hfmax)
+     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax)
       IMPLICIT NONE
 
@@ -60,5 +60,5 @@
 c    Argument added to account for subgrid gustiness :
 
-      REAL wmax(ngridmx), hfmax(ngridmx)
+      REAL wstar(ngridmx), hfmax(ngridmx), zi(ngridmx)
 
 c    Traceurs :
@@ -313,5 +313,5 @@
 c       ---------------------
 
-      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wmax,ptsrf,ph
+      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
      &             ,zcdv_true,zcdh_true)
 
@@ -319,14 +319,12 @@
 
         IF (callrichsl) THEN
-       zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(6.*hfmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
-       zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(10.*hfmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
-
-!       zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
-!       zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
-
-        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(6.*hfmax(:))**2)
+       zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
+     &  (14.75*sqrt(hfmax(:)) - 10.*hfmax(:) + 1.7*hfmax(:)**2)**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
+       zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
+     &  (20.6*hfmax(:) - 22.8*hfmax(:)**2 + 9.*hfmax(:)**4)**2)     ! LES comparisons. This parameter is linked to the thermals settings)
+
+        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
+     &  (14.75*sqrt(hfmax(:)) - 10.*hfmax(:) + 1.7*hfmax(:)**2)**2)
         tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
-!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)
-!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
 !        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
 !        tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
@@ -334,6 +332,4 @@
         zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance 
         zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
-!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
-!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
         ENDIF
 
