Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 495)
+++ /trunk/LMDZ.MARS/README	(revision 496)
@@ -1361,2 +1361,6 @@
 .true. will call functions geticecover to yield a nice evolving CO2 cap as measured by Titus (mostly useful for mesoscale)
 .false. is the default option, so this change is transparent to the casual user
+
+== 11/01/2012 == AC
+>> Corrections to the thermals scheme following latest revisions of the related paper. Replaced the surface layer interpolator by a more
+complete routine, that will ultimately become a post-processing utility (and disappear from libf) and a subroutine in the MCD.
Index: /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 496)
+++ /trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 496)
@@ -0,0 +1,443 @@
+      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)
+      IMPLICIT NONE
+!=======================================================================
+!
+!   Anlysis of the PBL from input temperature, wind field and thermals outputs.
+!
+!   -------  
+!
+!   Author: Arnaud Colaitis 09/01/12
+!   -------
+!
+!   Arguments:
+!   ----------
+!
+!   inputs:
+!   ------
+!     ngrid            size of the horizontal grid
+!     nlay             size of the vertical grid
+!     pz0(ngrid)       surface roughness length
+!     pg               gravity (m s -2)
+!     pz(ngrid,nlay)   height of layers
+!     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)
+!     hfmax(ngrid)     maximum vertical turbulent heat flux in thermals
+!     zmax(ngrid)      height reached by the thermals (pbl height)
+!     pts(ngrid)       surface temperature
+!     ph(ngrid,nlay)   potential temperature T*(p/ps)^kappa
+!     z_out            height of interpolation
+!
+!   outputs:
+!   ------
+!
+!     Teta_out(ngrid)  interpolated teta
+!     u_out(ngrid)     interpolated u
+!     ustar(ngrid)     friction velocity
+!     tstar(ngrid)     friction temperature
+!     wstar(ngrid)     free convection velocity
+!     L_mo(ngrid)      monin_obukhov length
+!
+!
+!=======================================================================
+!
+!-----------------------------------------------------------------------
+!   Declarations:
+!   -------------
+
+#include "comcstfi.h"
+#include "callkeys.h"
+
+!   Arguments:
+!   ----------
+
+      INTEGER, INTENT(IN) :: ngrid,nlay
+      REAL, INTENT(IN) :: pz0(ngrid)
+      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
+      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
+      REAL, INTENT(IN) :: wmax(ngrid),hfmax(ngrid),zmax(ngrid)
+      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
+      REAL, INTENT(IN) :: z_out
+      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)
+      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid),wstar(ngrid)
+      REAL, INTENT(OUT) :: L_mo(ngrid)
+
+!   Local:
+!   ------
+
+      INTEGER ig,k
+      REAL karman,nu
+      DATA karman,nu/.41,0.001/
+      SAVE karman,nu
+
+!    Local(2):
+!    ---------
+      
+      REAL zout
+      REAL rib(ngrid)  ! Bulk Richardson number
+      REAL fm(ngrid) ! stability function for momentum
+      REAL fh(ngrid) ! stability function for heat
+      REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
+          ! phim = 1+betam*zeta   or   (1-bm*zeta)**am
+          ! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
+      REAL betam, betah, alphah, bm, bh, lambda
+          ! ah and am are assumed to be -0.25 and -0.5 respectively
+      REAL cdn(ngrid),chn(ngrid)  ! neutral momentum and heat drag coefficient
+      REAL pz0t        ! initial thermal roughness length. (local)
+      REAL ric         ! critical richardson number
+      REAL reynolds(ngrid)    ! reynolds number for UBL
+      REAL prandtl(ngrid)     ! prandtl number for UBL
+      REAL pz0tcomp(ngrid)     ! computed z0t
+      REAL ite
+      REAL residual,zcd0,z1
+      REAL pcdv(ngrid),pcdh(ngrid)
+      REAL zu2(ngrid)                  ! Large-scale wind at first layer
+      REAL pbl_teta(ngrid)             ! mixed-layer potential temperature
+      INTEGER pbl_height_index(ngrid)  ! index of nearest vertical grid point for zmax
+      REAL dteta(ngrid,nlay),x(ngrid)  ! potential temperature gradient and z/zi
+      REAL dvhf(ngrid),dvvv(ngrid)     ! dimensionless vertical heat flux and 
+                                       ! dimensionless vertical velocity variance
+      REAL vhf(ngrid),vvv(ngrid)       ! vertical heat flux and vertical velocity variance
+      INTEGER ii(1)
+! temporary
+      INTEGER dimout
+
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+
+c Initialisation :
+
+      L_mo(:)=0.
+      ustar(:)=0.
+      tstar(:)=0.
+      zout=z_out
+      reynolds(:)=0.
+      pz0t = 0.
+      pz0tcomp(:) = 0.1*pz0(:)
+      rib(:)=0.
+      pcdv(:)=0.
+      pcdh(:)=0.
+
+      ! this formulation assumes alphah=1., implying betah=betam
+      ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
+
+      bm=16.            !UBL
+      bh=16.            !UBL
+      alphah=1.
+      betam=5.         !SBL
+      betah=5.         !SBL
+      lambda=(sqrt(bh/bm))/alphah
+      ric=betah/(betam**2)
+
+      DO ig=1,ngrid
+       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)
+
+       DO WHILE((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))
+
+         pz0t=pz0tcomp(ig)
+         IF (zu2(ig) .ne. 0.) THEN
+            ! 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)))
+     &        *(ph(ig,1)-pts(ig))/zu2(ig)
+         ELSE
+            print*,'warning, infinite Richardson at surface'
+            print*,pu(ig,1),pv(ig,1)
+            rib(ig) = ric
+         ENDIF
+
+         z1z0=pz(ig,1)/pz0(ig)
+         z1z0t=pz(ig,1)/pz0t
+
+         cdn(ig)=karman/log(z1z0)
+         cdn(ig)=cdn(ig)*cdn(ig)
+         chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) 
+
+         ! STABLE BOUNDARY LAYER :
+         IF (rib(ig) .gt. 0.) THEN
+            ! From D.E. England et al. (95)
+            prandtl(ig)=1.
+            if(rib(ig) .lt. ric) then
+               ! Assuming alphah=1. and bh=bm for stable conditions :
+               fm(ig)=((ric-rib(ig))/ric)**2
+               fh(ig)=((ric-rib(ig))/ric)**2
+            else
+               ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
+               fm(ig)=0.
+               fh(ig)=0.
+            endif
+         ! UNSTABLE BOUNDARY LAYER :
+         ELSE
+            ! From D.E. England et al. (95)
+            fm(ig)=sqrt(1.-lambda*bm*rib(ig))
+            fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*
+     &                     (1.-lambda*bm*rib(ig))**0.25
+            prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/
+     &             ((1.-lambda*bh*rib(ig))**0.5)
+         ENDIF
+ 
+        reynolds(ig)=karman*sqrt(fm(ig))
+     &              *sqrt(zu2(ig))
+     &              *pz0(ig)/(log(z1z0)*nu)
+        pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3*
+     &              (reynolds(ig)**0.25)*(prandtl(ig)**0.5))
+        residual = abs(pz0t-pz0tcomp(ig))
+        ite = ite+1
+ 
+       ENDDO  ! of while
+       pz0t=0.
+
+! Drag computation:
+
+         pcdv(ig)=cdn(ig)*fm(ig)
+         pcdh(ig)=chn(ig)*fh(ig)
+       
+      ENDDO ! of ngrid
+
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+
+
+! Large-scale wind at first layer (without gustiness) and
+! u* theta* computation
+
+      DO ig=1,ngrid
+         IF (rib(ig) .ge. ric) THEN
+           ustar(ig)=0.
+           tstar(ig)=0.
+         ELSE
+           ustar(ig)=sqrt(pcdv(ig))
+     &        *sqrt(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))
+           tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
+     &        /sqrt(pcdv(ig))
+         ENDIF
+      ENDDO
+
+! Monin Obukhov length:
+
+      DO ig=1,ngrid
+         IF (rib(ig) .gt. ric) THEN
+            L_mo(ig)=0.
+         ELSE
+            L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig))  !as defined here, L is positive for SBL, negative for UBL
+         ENDIF
+      ENDDO
+
+! Interpolation:
+
+      DO ig=1,ngrid
+        IF(zout .lt. pz0tcomp(ig)) THEN
+           u_out(ig)=0.
+           Teta_out(ig)=pts(ig)
+        ELSEIF (L_mo(ig) .gt. 0.) THEN
+           u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) + 
+     &            5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig))
+           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
+     &                            *log(zout/pz0tcomp(ig)) +
+     &            5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig)))
+     &                            *(zout-pz0tcomp(ig))
+        ELSEIF (L_mo(ig) .lt. 0.) THEN
+
+          IF(L_mo(ig) .gt. -1000.) THEN
+           
+           u_out(ig)=(ustar(ig)/karman)*(
+     &  2.*atan((1.-16.*zout/L_mo(ig))**0.25)
+     & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
+     & -2.*log(1.+(1.-16.*zout/L_mo(ig))**0.25)
+     & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
+     & -   log(1.+sqrt(1.-16.*zout/L_mo(ig)))
+     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
+     & +   log(zout/pz0(ig))
+     &                                  )
+
+           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
+     &  2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
+     & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig)))
+     & +   log(zout/pz0tcomp(ig))
+     &                                                        )
+
+          ELSE
+
+! We have to treat the case where L is very large and negative,
+! corresponding to a very nearly stable atmosphere (but not quite !)
+! i.e. teta* <0 and almost zero. We choose the cutoff
+! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
+! between the two expression for z/L = -1/1000 is -1.54324*10^-9
+! (we do that to avoid using r*4 precision, otherwise, we get -inf values)          
+
+           u_out(ig)=(ustar(ig)/karman)*(
+     &     (4./L_mo(ig))*(zout-pz0(ig))
+     &   + (20./(L_mo(ig))**2)*(zout**2-pz0(ig)**2)
+     &   + (160./(L_mo(ig))**3)*(zout**3-pz0(ig)**3)
+     &   + log(zout/pz0(ig))
+     &                                   )
+
+           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
+     &     (8./L_mo(ig))*(zout-pz0tcomp(ig))
+     &   + (48./(L_mo(ig))**2)*(zout**2-pz0tcomp(ig)**2)
+     &   + (1280./(3.*(L_mo(ig))**3))*(zout**3-pz0tcomp(ig)**3)
+     &   + log(zout/pz0tcomp(ig))
+     &                                                           )
+
+          ENDIF
+        ELSE
+           u_out(ig)=0.
+           Teta_out(ig)=pts(ig)
+        ENDIF
+        IF(zout .lt. pz0(ig)) THEN
+           u_out(ig)=0.
+        ENDIF
+      ENDDO
+
+! when using convective adjustment without thermals, a vertical potential temperature
+! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
+! interpolated at any height in the surface layer is theta at the first level.
+
+      IF ((.not.calltherm).and.(calladj)) THEN
+           Teta_out(:)=ph(:,1)
+      ENDIF
+
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+! PART III : WSTAR COMPUTATION
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+
+! Detection of the mixed-layer potential temperature
+! ------------
+
+! Nearest index for the pbl height
+
+      pbl_height_index(:)=1
+
+      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
+               pbl_height_index(ig)=k
+            ENDIF
+         ENDDO
+      ENDDO
+
+! Potential temperature gradient
+
+      dteta(:,nlay)=0.
+      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))
+         ENDDO
+      ENDDO
+
+! Computation of the pbl mixed layer temperature
+
+      DO ig=1, ngrid
+         ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig))))
+         pbl_teta(ig) = ph(ig,ii(1))
+      ENDDO
+
+! We follow Spiga et. al 2010 (QJRMS)
+! ------------
+
+      DO ig=1, ngrid
+         IF (zmax(ig) .gt. 0.) THEN
+            wstar(ig)=(pg*zmax(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.)
+         ELSE
+            wstar(ig)=0.
+         ENDIF
+      ENDDO
+
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+
+! We follow Spiga et. al 2010 (QJRMS)
+! ------------
+
+      DO ig=1, ngrid
+         IF (zmax(ig) .gt. 0.) THEN
+            x(ig) = zout/zmax(ig)
+         ELSE
+            x(ig) = 999.
+         ENDIF
+      ENDDO
+
+      DO ig=1, ngrid
+         ! dimensionless vertical heat flux
+         IF (x(ig) .le. 0.3) THEN
+            dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig)))
+     &                                       *exp(-4.61*x(ig))
+         ELSEIF (x(ig) .le. 1.) THEN
+            dvhf(ig) = -1.52*x(ig) + 1.24
+         ELSE
+            dvhf(ig) = 0.
+         ENDIF
+         ! dimensionless vertical velocity variance
+         IF (x(ig) .le. 1.) THEN
+            dvvv(ig) = 2.05*(x(ig)**(2./3.))*(1.-0.64*x(ig))**2
+         ELSE
+            dvvv(ig) = 0.
+         ENDIF
+      ENDDO
+
+      vhf(:) = dvhf(:)*hfmax(:)
+      vvv(:) = dvvv(:)*(wstar(:))**2
+
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+! OUTPUTS
+!------------------------------------------------------------------------
+!------------------------------------------------------------------------
+
+      IF (ngrid .eq. 1) THEN
+         dimout=0
+      ELSE
+         dimout=2
+      ENDIF
+
+         call WRITEDIAGFI(ngrid,'Teta_out',
+     &              'potential temperature at z_out','K',
+     &                         dimout,Teta_out)
+         call WRITEDIAGFI(ngrid,'u_out',
+     &              'horizontal velocity norm at z_out','m/s',
+     &                         dimout,u_out)
+         call WRITEDIAGFI(ngrid,'u_star',
+     &              'friction velocity','m/s',
+     &                         dimout,ustar)
+         call WRITEDIAGFI(ngrid,'teta_star',
+     &              'friction potential temperature','K',
+     &                         dimout,tstar)
+         call WRITEDIAGFI(ngrid,'L',
+     &              '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,'vvv',
+     &              'Vertical velocity variance at zout','m',
+     &                         dimout,vvv)
+         call WRITEDIAGFI(ngrid,'vhf',
+     &              'Vertical heat flux at zout','m',
+     &                         dimout,vhf)
+
+      RETURN
+      END
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 495)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 496)
@@ -333,5 +333,5 @@
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
       REAL, SAVE :: wmax_th(ngridmx)
-      REAL hfmax_th(ngridmx)
+      REAL, SAVE :: hfmax_th(ngridmx)
       REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
       REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
@@ -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)
+      REAL L_mo(ngridmx),wstar(ngridmx),vhf(ngridmx),vvv(ngridmx)
       REAL zu2(ngridmx)
 c=======================================================================
@@ -569,4 +569,6 @@
                  enddo
               endif
+           else
+              zdtnlte(:,:)=0.
            endif
 
@@ -746,5 +748,9 @@
              IF (zh(ig,1) .lt. tsurf(ig)) THEN
                wmax_th(ig)=1.
-             ENDIF        
+               hfmax_th(ig)=0.2
+             ELSE
+               wmax_th(ig)=0.
+               hfmax_th(ig)=0.
+             ENDIF      
           ENDDO
         ENDIF
@@ -762,5 +768,5 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
+     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh,hfmax_th)
 
 #ifdef MESOSCALE
@@ -853,4 +859,5 @@
         lmax_th(:)=0
         wmax_th(:)=0.
+        hfmax_th(:)=0.
         lmax_th_out(:)=0.
         end if
@@ -1803,27 +1810,9 @@
          z_out=0.
          if (calltherm .and. (z_out .gt. 0.)) then
-         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
-     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
-
-         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
-         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
-     &              'horizontal velocity norm','m/s',
-     &                         2,zu2)
-
-         call WRITEDIAGFI(ngridmx,'Teta_out',
-     &              'potential temperature at z_out','K',
-     &                         2,Teta_out)
-         call WRITEDIAGFI(ngridmx,'u_out',
-     &              'horizontal velocity norm at z_out','m/s',
-     &                         2,u_out)
-         call WRITEDIAGFI(ngridmx,'u*',
-     &              'friction velocity','m/s',
-     &                         2,ustar)
-         call WRITEDIAGFI(ngridmx,'teta*',
-     &              'friction potential temperature','K',
-     &                         2,tstar)
-         call WRITEDIAGFI(ngrid,'L',
-     &              'Monin Obukhov length','m',
-     &                         2,L_mo)
+
+         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)
+
          else
            if((.not. calltherm).and.(z_out .gt. 0.)) then
@@ -1914,24 +1903,9 @@
          z_out=0.
          if (calltherm .and. (z_out .gt. 0.)) then
-         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
-     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
-
-         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
-         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
-     &              'horizontal velocity norm','m/s',
-     &                         0,zu2)
-
-         call WRITEDIAGFI(ngridmx,'Teta_out',
-     &              'potential temperature at z_out','K',
-     &                         0,Teta_out)
-         call WRITEDIAGFI(ngridmx,'u_out',
-     &              'horizontal velocity norm at z_out','m/s',
-     &                         0,u_out)
-         call WRITEDIAGFI(ngridmx,'u*',
-     &              'friction velocity','m/s',
-     &                         0,ustar)
-         call WRITEDIAGFI(ngridmx,'teta*',
-     &              'friction potential temperature','K',
-     &                         0,tstar)
+
+         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)
+
          else
            if((.not. calltherm).and.(z_out .gt. 0.)) then
@@ -1947,4 +1921,7 @@
      &              'hauteur du thermique','point',
      &                         0,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'zmax_th',
+     &              'hauteur du thermique','m',
+     &                         0,zmax_th)
         call WRITEDIAGFI(ngridmx,'hfmax_th',
      &              'maximum TH heat flux','K.m/s',
Index: unk/LMDZ.MARS/libf/phymars/surflayer_interpol.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F	(revision 495)
+++ 	(revision )
@@ -1,371 +1,0 @@
-      SUBROUTINE surflayer_interpol(ngrid,nlay,pz0, 
-     & pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar,L_mo)
-      IMPLICIT NONE
-!=======================================================================
-!
-!   Subject: interpolation of temperature and velocity norm in the surface layer
-!   by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u*, teta*)
-!   -------  
-!
-!   Author: Arnaud Colaitis 05/08/11
-!   -------
-!
-!   Arguments:
-!   ----------
-!
-!   inputs:
-!   ------
-!     ngrid            size of the horizontal grid
-!     pg               gravity (m s -2)
-!     pz(ngrid,nlay)   height of layers
-!     pu(ngrid,nlay)   u component of the wind
-!     pv(ngrid,nlay)   v component of the wind
-!     pts(ngrid)       surface temperature
-!     ph(ngrid)        potential temperature T*(p/ps)^kappa
-!
-!
-!=======================================================================
-!
-!-----------------------------------------------------------------------
-!   Declarations:
-!   -------------
-
-#include "comcstfi.h"
-#include "callkeys.h"
-
-!   Arguments:
-!   ----------
-
-      INTEGER, INTENT(IN) :: ngrid,nlay
-      REAL, INTENT(IN) :: pz0(ngrid)
-      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
-      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
-      REAL, INTENT(IN) :: wmax(ngrid)
-      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
-      REAL, INTENT(IN) :: z_out                  ! output height (in m above surface)
-      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)! interpolated fields at z_out : potential temperature and norm(uv)
-      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta*
-      REAL, INTENT(OUT) :: L_mo(ngrid)                ! Monin-Obukhov length
-
-!   Local:
-!   ------
-
-      INTEGER ig
-
-      REAL karman,nu
-      DATA karman,nu/.41,0.001/
-      SAVE karman,nu
-
-!    Local(2):
-!    ---------
-      REAL zout
-
-      REAL rib(ngrid)  ! Bulk Richardson number
-      REAL fm(ngrid) ! stability function for momentum
-      REAL fh(ngrid) ! stability function for heat
-      REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
-
-! phim = 1+betam*zeta   or   (1-bm*zeta)**am
-! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
-      REAL betam, betah, alphah, bm, bh, lambda
-! ah and am are assumed to be -0.25 and -0.5 respectively
-
-      REAL cdn(ngrid),chn(ngrid)  ! neutral momentum and heat drag coefficient
-      REAL pz0t        ! initial thermal roughness length. (local)
-      REAL ric         ! critical richardson number
-      REAL reynolds(ngrid)    ! reynolds number for UBL
-      REAL prandtl(ngrid)     ! prandtl number for UBL
-      REAL pz0tcomp(ngrid)     ! computed z0t
-      REAL ite
-      REAL residual,zcd0,z1
-      REAL pcdv(ngrid),pcdh(ngrid)
-! For output :
-
-      REAL zu2(ngrid)                  ! Large-scale wind at first layer
-!-----------------------------------------------------------------------
-!   couche de surface:
-!   ------------------
-
-c Init :
-
-      L_mo(:)=0.
-      ustar(:)=0.
-      tstar(:)=0.
-      zout=z_out
-
-      reynolds(:)=0.
-
-! New formulation (AC) :
-
-! phim = 1+betam*zeta   or   (1-bm*zeta)**am
-! phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
-! am=-0.25, ah=-0.5
-
-      pz0t = 0.     ! for the sake of simplicity
-      pz0tcomp(:) = 0.1*pz0(:)
-      rib(:)=0.
-      pcdv(:)=0.
-      pcdh(:)=0.
-
-! this formulation assumes alphah=1., implying betah=betam
-! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
-      bm=16.            !UBL
-      bh=16.            !UBL
-      alphah=1.
-      betam=5.         !SBL
-      betah=5.         !SBL
-      lambda=(sqrt(bh/bm))/alphah
-      ric=betah/(betam**2)
-
-      DO ig=1,ngrid
-
-         ite=0.
-         residual=abs(pz0tcomp(ig)-pz0t)
-
-         do while((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))
-
-         pz0t=pz0tcomp(ig)
-
-         if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
-
-! Classical Richardson number formulation
-
-!         rib(ig) = (pg/ph(ig,1))*pz(ig,1)*(ph(ig,1)-pts(ig))
-!     &           /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
-
-! Richardson number formulation proposed by D.E. England et al. (1995)
-
-!        IF((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) .lt. 1.) 
-!     &            .and. (wmax(ig) .gt. 0.)) THEN
-          zu2(ig)=
-!     &  (MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
-     &  ( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
-!     &  pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
-!        ELSE
-!        zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
-!        ENDIF
-
-      if(.not.callrichsl) then
-         rib(ig) = (pg/ph(ig,1))
-     &      *sqrt(pz(ig,1)*pz0(ig))
-     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
-     &      *(ph(ig,1)-pts(ig))/(zu2(ig)+6.)
-      else
-          rib(ig) = (pg/ph(ig,1))
-!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
-     &      *sqrt(pz(ig,1)*pz0(ig))
-     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
-     &      *(ph(ig,1)-pts(ig))/(zu2(ig) + (0.5*wmax(ig))**2)
-
-!     &  /(MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
-!     &  /( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
-
-      endif
-
-         else
-         print*,'warning, infinite Richardson at surface'
-         print*,pu(ig,1),pv(ig,1)
-         rib(ig) = ric          ! traiter ce cas !
-         endif
-
-         z1z0=pz(ig,1)/pz0(ig)
-         z1z0t=pz(ig,1)/pz0t
-
-         cdn(ig)=karman/log(z1z0)
-         cdn(ig)=cdn(ig)*cdn(ig)
-         chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) 
-
-! Stable case :
-      if (rib(ig) .gt. 0.) then
-! From D.E. England et al. (95)
-      prandtl(ig)=1.
-         if(rib(ig) .lt. ric) then
-! Assuming alphah=1. and bh=bm for stable conditions :
-            fm(ig)=((ric-rib(ig))/ric)**2
-            fh(ig)=((ric-rib(ig))/ric)**2
-         else
-! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
-            fm(ig)=0.
-            fh(ig)=0.
-         endif
-! Unstable case :
-      else
-! From D.E. England et al. (95)
-         fm(ig)=sqrt(1.-lambda*bm*rib(ig))
-         fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*
-     &                     (1.-lambda*bm*rib(ig))**0.25
-         prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/
-     &             ((1.-lambda*bh*rib(ig))**0.5)
-      endif
-
-       reynolds(ig)=karman*sqrt(fm(ig))
-     &              *sqrt(zu2(ig)+(0.5*wmax(ig))**2)
-!    &               *sqrt(pu(ig,1)**2 + pv(ig,1)**2)
-     &              *pz0(ig)/(log(z1z0)*nu)
-       pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3*
-     &              (reynolds(ig)**0.25)*(prandtl(ig)**0.5))
-
-          
-      residual = abs(pz0t-pz0tcomp(ig))
-      ite = ite+1
-!      if(ite .eq. 10) then
-!        print*, 'iteration max reached'
-!      endif
-!      print*, "iteration nnumber, residual",ite,residual
-
-      enddo  ! of while
-
-          pz0t=0.
-
-! Drag computation :
-
-         pcdv(ig)=cdn(ig)*fm(ig)
-         pcdh(ig)=chn(ig)*fh(ig)
-       
-      ENDDO
-
-!      endif !of if callrichsl
-
-! Large-scale wind at first layer (without gustiness) and
-! u* theta* computation
-      DO ig=1,ngrid
-
-         if (rib(ig) .ge. ric) then
-           ustar(ig)=0.
-           tstar(ig)=0.
-         else
-
-!           ustar(ig)=karman*sqrt(fm(ig)*zu2(ig))/(log(pz(ig,1)/pz0(ig)))
-!           tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
-!     &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
-
-!simpler definition of u* and teta*, equivalent to the formula above :
-
-            ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))
-            tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
-     &     /sqrt(pcdv(ig))
-
-           if ((tstar(ig) .lt. -50) .and. callrichsl) then
-             print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
-     &               ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
-           endif 
-         endif
-      ENDDO
-
-! Monin Obukhov length :
-
-      DO ig=1,ngrid
-         if (rib(ig) .gt. ric) then
-            L_mo(ig)=0.
-         else
-            L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig))  !as defined here, L is positive for SBL, negative for UBL
-         endif
-      ENDDO
-
-      DO ig=1,ngrid
-      IF(zout .ge. pz(ig,1)) THEN
-           zout=1.
-           print*, 'Warning, z_out should be less than the first 
-     & vertical grid-point'
-           print*, 'z1 =',pz(ig,1)
-           print*, 'z_out=',z_out
-           print*, 'z_out has been changed to 1m 
-     & and computation is resumed'
-      ENDIF
-
-      IF(zout .lt. pz0tcomp(ig)) THEN
-           zout=pz0tcomp(ig)
-           print*, 'Warning, z_out should be more than the thermal 
-     & roughness length'
-           print*, 'z0 =',pz0tcomp(ig)
-           print*, 'z_out=',z_out
-           print*, 'z_out has been changed to z0t
-     & and computation is resumed'
-      ENDIF
-      ENDDO
-
-      print*, 'interpolation of u and teta at z_out=',zout
-
-      DO ig=1,ngrid
-        IF (L_mo(ig) .gt. 0.) THEN
-           u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) + 
-     &            5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig))
-           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
-     &                            *log(zout/pz0tcomp(ig)) +
-     &            5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig)))
-     &                            *(zout-pz0tcomp(ig))
-        ELSEIF (L_mo(ig) .lt. 0.) THEN
-
-          IF(L_mo(ig) .gt. -1000.) THEN
-           
-           u_out(ig)=(ustar(ig)/karman)*(
-     &  2.*atan((1.-16.*zout/L_mo(ig))**0.25)
-     & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
-     & -2.*log(1.+(1.-16.*zout/L_mo(ig))**0.25)
-     & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
-     & -   log(1.+sqrt(1.-16.*zout/L_mo(ig)))
-     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
-     & +   log(zout/pz0(ig))
-     &                                  )
-
-           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
-     &  2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
-     & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig)))
-     & +   log(zout/pz0tcomp(ig))
-     &                                                        )
-
-          ELSE
-
-! We have to treat the case where L is very large and negative,
-! corresponding to a very nearly stable atmosphere (but not quite !)
-! i.e. teta* <0 and almost zero. We choose the cutoff
-! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
-! between the two expression for z/L = -1/1000 is -1.54324*10^-9
-! (we do that to avoid using r*4 precision, otherwise, we get -inf values)          
-
-           u_out(ig)=(ustar(ig)/karman)*(
-     &     (4./L_mo(ig))*(zout-pz0(ig))
-     &   + (20./(L_mo(ig))**2)*(zout**2-pz0(ig)**2)
-     &   + (160./(L_mo(ig))**3)*(zout**3-pz0(ig)**3)
-     &   + log(zout/pz0(ig))
-     &                                   )
-
-           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
-     &     (8./L_mo(ig))*(zout-pz0tcomp(ig))
-     &   + (48./(L_mo(ig))**2)*(zout**2-pz0tcomp(ig)**2)
-     &   + (1280./(3.*(L_mo(ig))**3))*(zout**3-pz0tcomp(ig)**3)
-     &   + log(zout/pz0tcomp(ig))
-     &                                                           )
-
-          ENDIF
-        ELSE
-           u_out(ig)=0.
-           Teta_out(ig)=pts(ig)
-        ENDIF
-      ENDDO
-
-! when using convective adjustment without thermals, a vertical potential temperature
-! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
-! interpolated at any height in the surface layer is theta at the first level.
-
-      IF ((.not.calltherm).and.(calladj)) THEN
-           Teta_out(:)=ph(:,1)
-      ENDIF
-
-! Usefull diagnostics for the interpolation routine :
-
-!
-!         call WRITEDIAGFI(ngrid,'Ri',
-!     &              'Richardson','m',
-!     &                         2,rib)
-!
-!         call WRITEDIAGFI(ngrid,'z0T',
-!     &              'thermal roughness length','m',
-!     &                         2,pz0tcomp)
-!         call WRITEDIAGFI(ngrid,'z0',
-!     &              'roughness length','m',
-!     &                         2,pz0)
-
-
-      RETURN
-      END
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 495)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 496)
@@ -111,5 +111,5 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL active(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd
+      REAL a1,b1,ae,be,ad,bd,fdfu
       INTEGER tic
 
@@ -304,8 +304,23 @@
       linter(:)=1.
 
+! --------------------------------------------------------------------------
+! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------
+! --------------  see thermiques.pro and getfit.py -------------------------
+
 !      a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6  ! svn baseline
 
-      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57 !improved fits
-      ad = 0.0005114  ; bd = -0.662
+! Using broad downdraft selection
+!      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57
+!      ad = 0.0005114  ; bd = -0.662
+!      fdfu = -1.9
+
+! Using conditional sampling downdraft selection
+      a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421
+      ad = 0.00048088  ; bd = -0.6697
+      fdfu = -1.3
+
+! --------------------------------------------------------------------------
+! --------------------------------------------------------------------------
+! --------------------------------------------------------------------------
 
 ! Initialisation des variables entieres
@@ -513,5 +528,5 @@
 !---------------------------------------------------------------------------
 
-      DO tic=0,6  ! internal convergence loop
+      DO tic=0,1  ! internal convergence loop
       activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -640,5 +655,4 @@
 ! ================= FIN PLUME ===============================================
 ! ===========================================================================
-
 
 ! ===========================================================================
@@ -1089,5 +1103,6 @@
          do l=1,lmax(ig)
               if(zlay(ig,l) .le. 0.7*zmax(ig)) then
-              fm_down(ig,l) =-1.9*fm(ig,l)
+              fm_down(ig,l) =fm(ig,l)* &
+     &      max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
               endif
 
@@ -1102,14 +1117,4 @@
             endif
 
-!             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig))/19.231 + 0.9938)
-!             elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig)-0.075)/187.931 + 0.9982)
-!             else
-!!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig) -0.60)/(-1333) + 1.00025)
-!!             else
-!!          ztvd(ig,l)=ztv(ig,l)
-!            endif
 
          enddo
@@ -1126,5 +1131,6 @@
          if(lmax(ig) .gt. 1) then
 ! No downdraft in the very-near surface layer, we begin at k=3
-         do k=2,lmax(ig)
+ 
+         do k=3,lmax(ig)
             zdthladj_down(ig,k)=(1./masse(ig,k))*(fm_down(ig,k+1)*ztv(ig,k+1)- &
      & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
@@ -1332,4 +1338,5 @@
          do ig=1,ngridmx
            pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
+!           pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)
          enddo
       enddo
Index: /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 495)
+++ /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 496)
@@ -148,6 +148,6 @@
 !         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wmax(ig)**2)
 !         zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
-!         zu2=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
+         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
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 495)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 496)
@@ -5,5 +5,5 @@
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
-     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
+     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true,hfmax)
       IMPLICIT NONE
 
@@ -44,5 +44,5 @@
       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
       REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
-      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay),pt(ngrid,nlay)
       REAL ptsrf(ngrid),pemis(ngrid)
       REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
@@ -60,5 +60,5 @@
 c    Argument added to account for subgrid gustiness :
 
-      REAL wmax(ngridmx)
+      REAL wmax(ngridmx), hfmax(ngridmx)
 
 c    Traceurs :
@@ -78,4 +78,5 @@
       REAL z4st,zdplanck(ngridmx)
       REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
+      REAL zkq(ngridmx,nlayermx+1)
       REAL zcdv(ngridmx),zcdh(ngridmx)
       REAL zcdv_true(ngridmx),zcdh_true(ngridmx)    ! drag coeff are used by the LES to recompute u* and hfx
@@ -114,5 +115,4 @@
       REAL kmixmin
 
-
 c    Mass-variation scheme :
 c    ~~~~~~~
@@ -124,7 +124,6 @@
       REAL pdtc(ngrid,nlayermx)
       REAL zhs(ngridmx,nlayermx)
-      REAL cpice,ccond
-      SAVE ccond,cpice
-      DATA cpice /1000./
+      REAL ccond
+      SAVE ccond
 
 c     Theta_m formulation for mass-variation scheme :
@@ -320,6 +319,12 @@
 
         IF (callrichsl) THEN
-        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)
+       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)
+        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
 !        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)
 !        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
@@ -342,7 +347,7 @@
 !     &              'heat drag','no units',
 !     &                         2,zcdh_true)
-!        call WRITEDIAGFI(ngridmx,'u*',
+!        call WRITEDIAGFI(ngridmx,'ust',
 !     &              'friction velocity','m/s',2,ust)
-!       call WRITEDIAGFI(ngridmx,'T*',
+!       call WRITEDIAGFI(ngridmx,'tst',
 !     &              'friction temperature','K',2,tst)
 !        call WRITEDIAGFI(ngridmx,'rm-1',
@@ -355,8 +360,13 @@
 c    ** schema de diffusion turbulente dans la couche limite
 c       ---------------------------------------------------- 
-
-        CALL vdif_kc(ptimestep,g,pzlev,pzlay
+!
+       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
+!     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
+!     s   ,8)
 
       if ((doubleq).and.(ngrid.eq.1)) then
