Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 289)
+++ trunk/LMDZ.MARS/README	(revision 290)
@@ -965,2 +965,43 @@
 with the new mixing layer scheme (thermals). And anyway it is a much
 better scheme than the previous one.
+
+== 09/09/11 == AC
+
+ This is a major update for thermals and richardson layer parametrization. This update stabilizes thermals (further
+ studies might show that we can reduce the value of nsplit in gcm. To be tested.), and prevent downdrafts from descending into
+ the surface layer, which was acting as a cold feedback on the thermals. The Richardson surface layer now features
+ different gustiness coefficients for Richardson, heat and momentum so that u* and t* are correctly represented.
+ 
+ Upcoming updates will change surflayer_interpol.F90 to implement those changes in the interpolation scheme as well.
+ 
+ ***************************
+ IMPORTANT : several parts of the vdifc code might want to use these new definitions for gustiness, u* and t*. exemple : dust devil routines
+ that recompute u* ? lifting routines ? TODO !
+ **************************
+ 
+ 
+ M             289   libf/phymars/thermcell_main_mars.F90
+ ^-----------------> Added iterations to the velocity / buoyancy / entrainment / detrainment
+                     computation to ensure correct convergence. Iteration number is for now set to
+                     4, which is probably too much. This will be changed once several cases are tested.
+                     The minimum is probably 2 iterations.
+ 
+ M             289   libf/phymars/vdifc.F
+ ^-----------------> Subgrid gustiness parametrization now uses different gustiness beta coefficients
+                     for heat and momentum. Comparisons with LES at Exomars landing site, matching u*
+                     and teta* suggests values of beta=0.7 for momentum and beta=1.2 for heat, where
+                     the formula for large scale horizontal winds in the first layer is :
+                          U0^2 = pu^2 + pv^2 + (beta*wmax_th)^2
+ 
+ M             289   libf/phymars/vdif_cd.F
+ ^-----------------> LES data suggests that Richardson number distribution during daytime has a very large
+                     standart deviation due to highly negative Richardsons on several gridpoints. As a consequence,
+                     the mean Richardson in daytime in the LES is much more negative than in LES. In the gcm,
+                     parametrized subgrid turbulence prevents such cases, which can be dangerous in nearly unstable conditions.
+                     Hence, we use beta=0.5 instead of one, so that we keep the anti-crazy-hfx function of beta and we increase the
+                     norm of negative Richardsons in general for daytime conditions. This is set in conjonction with
+                     beta settings for heat and momentum in vdifc.
+ 
+ M             289   libf/phymars/meso_inc/meso_inc_les.F
+ ^-----------------> HFX and USTM computations now uses the different betas for heat and momentum.
+
Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F	(revision 289)
+++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F	(revision 290)
@@ -20,11 +20,10 @@
             hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp
      &        *sqrt(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
-     &        + (1.0*wmax_th(ig))**2)
+     &        + (1.2*wmax_th(ig))**2)
      &        *zcdh(ig)*(tsurf(ig)-zh(ig,1))
-
 
 ! New SL parametrization, ust is more accurately computed in vdif_cd :
             ust(ig) = sqrt(zcdv(ig)*
-     &   (pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (1.0*wmax_th(ig))**2)
+     &   (pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.7*wmax_th(ig))**2)
      &                     )
 
Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 289)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 290)
@@ -114,4 +114,5 @@
       LOGICAL active(ngridmx),activetmp(ngridmx)
       REAL a1,b1,ae,be,ad,bd
+      INTEGER tic
 
 ! ==========================================
@@ -334,4 +335,5 @@
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
      &                       *sqrt(zlev(ig,2))
+!     &                       /sqrt(zlev(ig,2))
 !      &                       *zlev(ig,2)
                lalim(ig)=2
@@ -381,8 +383,10 @@
 ! dans un panache conservatif
           f_star(ig,1)=0.
+          
           f_star(ig,2)=alim_star(ig,1)
+
           zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
       &                     *(zlev(ig,2)-zlev(ig,1))  &
-      &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
+      &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
           w_est(ig,2)=zw2(ig,2)
 
@@ -420,5 +424,5 @@
 
 !                if(l .lt. lalim(ig)) then
-!                ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
+!               ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
 !     &            alim_star(ig,l)*ztv(ig,l))  &
 !     &            /(f_star(ig,l)+alim_star(ig,l))
@@ -522,4 +526,6 @@
 !---------------------------------------------------------------------------
 
+      DO tic=0,3
+
       activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -547,4 +553,52 @@
       enddo
 
+! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
+
+
+      do ig=1,ngridmx
+        if (activetmp(ig)) then
+
+          zw2m=zw2(ig,l+1)
+          if(zw2m .gt. 0) then
+          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
+          entr_star(ig,l)=f_star(ig,l)*zdz*  &
+        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+          else
+          entr_star(ig,l)=0.
+          endif
+
+          if(zbuoy(ig,l) .gt. 0.) then
+             if(l .lt. lalim(ig)) then
+                detr_star(ig,l)=0.
+             else
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  ad
+             endif
+          else
+          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
+                &     bd*zbuoy(ig,l)/zw2m
+          endif
+          else
+          entr_star(ig,l)=0.
+          detr_star(ig,l)=0.
+          endif
+
+! En dessous de lalim, on prend le max de alim_star et entr_star pour
+! alim_star et 0 sinon
+
+        if (l.lt.lalim(ig)) then
+          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
+          entr_star(ig,l)=0.
+        endif
+
+! Calcul du flux montant normalise
+
+      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
+     &              -detr_star(ig,l)
+
+      endif
+      enddo
+
+      ENDDO
 !---------------------------------------------------------------------------
 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
@@ -1121,5 +1175,6 @@
       do ig=1,ngridmx
          if(lmax(ig) .gt. 1) then
-         do k=1,lmax(ig)
+! No downdraft in the very-near surface layer, we begin at k=3
+         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))
Index: trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 289)
+++ trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 290)
@@ -149,5 +149,5 @@
 !         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) + (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 289)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 290)
@@ -234,16 +234,14 @@
      &             ,zcdv_true,zcdh_true)
 
-
-        IF (callrichsl .eq. .true.) then
-           zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
-     &             +(1.*wmax(:))**2        ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
-                                            ! LES comparisons. This parameter is linked to the thermals settings)
+        zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
+
+        IF (callrichsl .eq. .true.) 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)
         ELSE
-           zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
-        ENDIF        
-
         zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance 
         zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
-                                             ! these are the quantities to be looked at when comparing surface layers of different models
+        ENDIF
+
 
 ! Some usefull diagnostics for the new surface layer parametrization :
@@ -257,8 +255,10 @@
 !        call WRITEDIAGFI(ngridmx,'u*',
 !     &              'friction velocity','m/s',
-!     &                         2,sqrt(zcdv_true(:))*sqrt(zu2(:)))
+!     &             2,sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2))
 !        call WRITEDIAGFI(ngridmx,'T*',
 !     &              'friction temperature','K',
-!     &          2,-zcdh_true(:)*(ptsrf(:)-ph(:,1))/sqrt(zcdv_true(:)))
+!     &   2,zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)*
+!     &   -(ptsrf(:)-ph(:,1))
+!     & /(sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)))
 !        call WRITEDIAGFI(ngridmx,'rm-1',
 !     &              'aerodyn momentum conductance','m/s',
