Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 511)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 512)
@@ -59,5 +59,5 @@
       REAL detr_therm(ngridmx,nlayermx)
       REAL zw2(ngridmx,nlayermx+1)
-      REAL fraca(ngridmx,nlayermx+1)
+      REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
       REAL ztla(ngridmx,nlayermx)
       REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
@@ -124,4 +124,5 @@
       zw2(:,:)=0.
       fraca(:,:)=0.
+      zfraca(:,:)=0.
       if (tracer) then
          pdq_th(:,:,:)=0.
@@ -264,4 +265,5 @@
             detr_therm(:,:)=detr_therm(:,:)  &
      &       +zdetr_therm(:,:)*fact
+            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
 
             heatFlux(:,:)=heatFlux(:,:) &
@@ -273,4 +275,5 @@
             buoyancyEst(:,:)=buoyancyEst(:,:) &
      &       +zbuoyancyEst(:,:)*fact
+  
 
             zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 511)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 512)
@@ -113,5 +113,5 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL activecell(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd,fdfu
+      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv
       INTEGER tic
 
@@ -379,4 +379,10 @@
       fdfu = -1.3
 
+! Best configuration for 222 levels:
+      b1=0.
+      a1=0.7*a1
+      a1inv=0.25*a1
+      b1inv=0.0002
+
 ! --------------------------------------------------------------------------
 ! --------------------------------------------------------------------------
@@ -406,5 +412,5 @@
 !        do l=2,4
          do ig=1,ngridmx
-           if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
+           if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.) .and. (zlev(ig,l+1) .lt. 1000.)) then
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
      &                       *sqrt(zlev(ig,l+1))
@@ -499,5 +505,5 @@
      & -2.*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be)
                 else
-                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1)
+                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1inv)
                 endif
                 if (w_est(ig,l+1).lt.0.) then
@@ -518,6 +524,6 @@
           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)
-         &  MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
+        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -540,6 +546,6 @@
 ! new param from continuity eq with a fit on dfdz
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
-!            &  ad
-             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+            &  ad
+!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
 
 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
@@ -553,6 +559,7 @@
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
-!                &     bd*zbuoy(ig,l)/zw2m
-             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+                &     bd*zbuoy(ig,l)/zw2m
+!             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+!             &  Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
 
 
@@ -573,8 +580,8 @@
 ! alim_star et 0 sinon
 
-        if (l.lt.lalim(ig)) then
+       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
+          entr_star(ig,l)=0. 
+       endif
 
 ! Calcul du flux montant normalise
@@ -591,5 +598,5 @@
 !---------------------------------------------------------------------------
 
-      DO tic=0,1  ! internal convergence loop
+      DO tic=0,0  ! internal convergence loop
       activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -614,5 +621,5 @@
      & 2.*zdz*zw2(ig,l)*b1-2.*zdz*zw2(ig,l)*ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be)
            else
-           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1)
+           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1inv)
            endif
       endif
@@ -630,4 +637,5 @@
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
         &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -640,8 +648,12 @@
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
             &  ad
+!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
+
              endif
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
                 &     bd*zbuoy(ig,l)/zw2m
+!             &  Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
+
           endif
           else
@@ -736,5 +748,5 @@
          do l=nlayermx,lalim(ig)+1,-1
             if (zw2(ig,l).le.1.e-10) then
-               lmax(ig)=l-1
+               lmax(ig)=l-1 
             endif
          enddo
@@ -1165,19 +1177,43 @@
          if (lmax(ig) .gt. 1) then
          do l=1,lmax(ig)
-              if(zlay(ig,l) .le. 0.7*zmax(ig)) then
+!              if(zlay(ig,l) .le. 0.8*zmax(ig)) then
+              if(zlay(ig,l) .le. zmax(ig)) then
               fm_down(ig,l) =fm(ig,l)* &
      &      max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
               endif
 
-             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
-             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
-             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
+!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+!            if(zlay(ig,l) .le. 0.6*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
+!            elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
+!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+
+!            if(zlay(ig,l) .le. 0.8*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
+!            elseif(zlay(ig,l) .le. zmax(ig)) then
+!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+
+            if(zlay(ig,l) .le. zmax(ig)) then
+          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
              else
           ztvd(ig,l)=ztv(ig,l)
             endif
-
 
          enddo
@@ -1432,4 +1468,5 @@
          do ig=1,ngridmx
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
+!         pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l)
          enddo
       enddo
