Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 1130)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 1212)
@@ -37,5 +37,5 @@
      &                  ,pu,pv,pt,pq,pq2  &
      &                  ,pdtadj,pdqadj  &
-     &                  ,fm,entr,detr,lmax,zmax  &
+     &                  ,fm,entr,detr,lmax,zmax,limz  &
      &                  ,zw2,fraca &
      &                  ,zpopsk,heatFlux,heatFlux_down &
@@ -43,5 +43,7 @@
 
       USE comtherm_h
+#ifndef MESOSCALE
       use planetwide_mod, only: planetwide_maxval
+#endif
 
       IMPLICIT NONE
@@ -85,4 +87,6 @@
       REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlayer) ! interface heat flux from downdraft
 
+      INTEGER, INTENT(OUT) :: limz ! limit vertical index for integration
+
 ! ============== LOCAL ================
       REAL :: pdqadj(ngrid,nlayer,nq) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop)
@@ -153,5 +157,4 @@
       REAL denom(ngrid)
       REAL zlevinter(ngrid)
-      INTEGER zlmax
 
 ! =========================================
@@ -765,12 +768,16 @@
 ! ===========================================================================
 
-      !zlmax=MAXVAL(lmax(:))+2 ! OK, but in serial mode only; use planet
-      call planetwide_maxval(lmax,zlmax)
-      zlmax=zlmax+2
+#ifdef MESOSCALE
+      limz= nlayer-5 ! the most important is limz > max(PBLheight)+2
+                     ! nlayer-5 is more than enough!
+#else
+      call planetwide_maxval(lmax,limz)
+      limz=limz+2
+#endif
       
-      if (zlmax .ge. nlayer) then
+      if (limz .ge. nlayer) then
         print*,'thermals have reached last layer of the model'
         print*,'this is not good !'
-        zlmax=nlayer
+        limz=nlayer
       endif
 ! alim_star_clos is the source profile used for closure. It consists of the
@@ -794,9 +801,13 @@
 
 ! llmax is the index of the heighest thermal in the simulation domain
-      !llmax=1
-      !do ig=1,ngrid
-      !   if (lalim(ig)>llmax) llmax=lalim(ig)
-      !enddo
+#ifdef MESOSCALE
+      !! AS: THIS IS PARALLEL SENSITIVE!!!!! to be corrected?
+      llmax=1
+      do ig=1,ngrid
+         if (lalim(ig)>llmax) llmax=lalim(ig)
+      enddo
+#else
       call planetwide_maxval(lalim,llmax)
+#endif
 
 ! Integral of a**2/(rho* Delta z), see equation 13 of appendix 4.2 in paper
@@ -856,5 +867,5 @@
 !-------------------------------------------------------------------------
 
-      do l=1,zlmax
+      do l=1,limz
          entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
          detr(:,l)=f(:)*detr_star(:,l)
@@ -863,5 +874,5 @@
 ! Reconstruct the updraft mass flux everywhere
 
-      do l=1,zlmax
+      do l=1,limz
          do ig=1,ngrid
             if (l.lt.lmax(ig)) then
@@ -890,5 +901,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      do l=1,zlmax !loop on the levels
+      do l=1,limz !loop on the levels
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -992,5 +1003,5 @@
 !-----------------------------------------------------------------------
 
-      do l=1,zlmax
+      do l=1,limz
          do ig=1,ngrid
             eee0=entr(ig,l)
@@ -1140,5 +1151,5 @@
 !------------------------------------------------------------------
       fraca(:,:)=0.
-      do l=2,zlmax
+      do l=2,limz
          do ig=1,ngrid
             if (zw2(ig,l).gt.1.e-10) then
@@ -1161,5 +1172,5 @@
       if (igcm_co2.ne.0) then
       detrmod(:,:)=0.
-      do k=1,zlmax
+      do k=1,limz
          do ig=1,ngrid
             detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
@@ -1174,9 +1185,9 @@
       call thermcell_dqup(ngrid,nlayer,ptimestep     &
      &     ,fm,entr,detrmod,  &
-     &    masse,pq(:,:,igcm_co2),pdqadj(:,:,igcm_co2),ndt,zlmax)
+     &    masse,pq(:,:,igcm_co2),pdqadj(:,:,igcm_co2),ndt,limz)
 
 ! Compute the ratio between theta and theta_m
       
-       do l=1,zlmax
+       do l=1,limz
           do ig=1,ngrid
              ratiom(ig,l)=1./(A*(pq(ig,l,igcm_co2)+pdqadj(ig,l,igcm_co2)*ptimestep)+B)
@@ -1191,5 +1202,5 @@
 
       pdtadj(:,:)=0.
-      do l=1,zlmax
+      do l=1,limz
          do ig=1,ngrid
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
@@ -1225,5 +1236,5 @@
         buoyancyEst(:,:)=0.
         heatFlux_down(:,:)=0.
-      do l=1,zlmax
+      do l=1,limz
        do ig=1,ngrid
         heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
