Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 184)
+++ /trunk/LMDZ.MARS/README	(revision 185)
@@ -763,2 +763,9 @@
 - correct bug (introduced previously) in lect_start_archive.F on loop
   boundaries for soil temperature.
+
+== 17/06/2011 == AC
+- Added new settings for the Martian thermals from new LES observations
+- Revamped thermcell's module variables to allow it's removal
+- Minor changes in physiq and meso_physiq for the call to thermals
+- Switched from dynamic to static memory allocation for all thermals variable 
+to gain computation speed
Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 185)
@@ -2,9 +2,9 @@
 ! AC 2011-01-05
 !
-      SUBROUTINE calltherm_interface (ngrid,nlayer,firstcall, &
+      SUBROUTINE calltherm_interface (firstcall, &
      & long,lati,zzlev,zzlay, &
      & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
-     & pplay,pplev,pphi,nq,zpopsk, &
-     & pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,pbl_dtke,hfmax,wmax)
+     & pplay,pplev,pphi,zpopsk, &
+     & pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,pbl_dtke,hfmax,wmax)
 
        USE ioipsl_getincom
@@ -12,22 +12,24 @@
       implicit none
 #include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
 !--------------------------------------------------------
 ! Variables d'entree
 !--------------------------------------------------------
 
-      INTEGER, INTENT(IN) :: ngrid,nlayer,nq
       REAL, INTENT(IN) :: ptimestep
-      REAL, INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL, INTENT(IN) :: pphi(ngrid,nlayer)
-      REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer)
-      REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
-      REAL, INTENT(IN) :: zzlay(ngrid,nlayer)
-      REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) 
+      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)
+      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1) 
       LOGICAL, INTENT(IN) :: firstcall
-      REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer)
-      REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq),pdt(ngrid,nlayer)
-      REAL, INTENT(IN) :: q2(ngrid,nlayer+1)
-      REAL, INTENT(IN) :: long(ngrid),lati(ngrid)
-      REAL, INTENT(IN) :: zpopsk(ngrid,nlayer)
+      REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
+      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
+      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
 
 !--------------------------------------------------------
@@ -35,25 +37,26 @@
 !--------------------------------------------------------
 
-      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
-      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
-      INTEGER lmax_th(ngrid)
-      REAL pbl_dtke(ngrid,nlayer+1)
+      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
+      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
+      INTEGER lmax_th(ngridmx)
+      REAL zmax_th(ngridmx)
+      REAL pbl_dtke(ngridmx,nlayermx+1)
 
 !--------------------------------------------------------
 ! Variables du thermique
 !--------------------------------------------------------
-      REAL u_seri(ngrid,nlayer), v_seri(ngrid,nlayer)
-      REAL t_seri(ngrid,nlayer)
-      REAL d_t_ajs(ngrid,nlayer)
-      REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
-      REAL d_v_ajs(ngrid,nlayer) 
-      REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer)
-      REAL detr_therm(ngrid,nlayer)
-      REAL zw2(ngrid,nlayer+1)
-      REAL fraca(ngrid,nlayer+1)
-      REAL ztla(ngrid,nlayer)
-      REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq)
-      REAL dq_therm(ngrid,nlayer), dq_thermdown(ngrid,nlayer)
-      REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer)
+      REAL u_seri(ngridmx,nlayermx), v_seri(ngridmx,nlayermx)
+      REAL t_seri(ngridmx,nlayermx)
+      REAL d_t_ajs(ngridmx,nlayermx)
+      REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
+      REAL d_v_ajs(ngridmx,nlayermx) 
+      REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
+      REAL detr_therm(ngridmx,nlayermx)
+      REAL zw2(ngridmx,nlayermx+1)
+      REAL fraca(ngridmx,nlayermx+1)
+      REAL ztla(ngridmx,nlayermx)
+      REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
+      REAL dq_therm(ngridmx,nlayermx), dq_thermdown(ngridmx,nlayermx)
+      REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
 
       LOGICAL qtransport_thermals,dtke_thermals
@@ -63,9 +66,11 @@
 ! Variable de diagnostique : flux de chaleur vertical
 
-      REAL heatFlux(ngrid,nlayer)
-      REAL heatFlux_down(ngrid,nlayer)
-      REAL buoyancyOut(ngrid,nlayer)
-      REAL buoyancyEst(ngrid,nlayer)
-      REAL hfmax(ngrid),wmax(ngrid)
+      REAL heatFlux(ngridmx,nlayermx)
+      REAL heatFlux_down(ngridmx,nlayermx)
+      REAL buoyancyOut(ngridmx,nlayermx)
+      REAL buoyancyEst(ngridmx,nlayermx)
+      REAL hfmax(ngridmx),wmax(ngridmx)
+
+      REAL tstart,tstop
 
 !---------------------------------------------------------
@@ -125,18 +130,22 @@
          if(dtke_thermals) then
 
-         DO l=1,nlayer
+         DO l=1,nlayermx
               q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
          ENDDO
          endif
 
-         CALL calltherm_mars(ngrid,nlayer,ptimestep,nq,zzlev,zzlay &
+         call cpu_time(tstart)
+
+         CALL calltherm_mars(ptimestep,zzlev,zzlay &
      &      ,pplay,pplev,pphi &
      &      ,u_seri,v_seri,t_seri,pq_therm, q2_therm &
      &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, dq2_therm &
      &      ,fm_therm,entr_therm,detr_therm &
-     &      ,lmax_th &
+     &      ,lmax_th,zmax_th &
      &      ,zw2,fraca &
      &      ,zpopsk,ztla,heatFlux,heatFlux_down &
      &      ,buoyancyOut,buoyancyEst,hfmax,wmax)
+          call cpu_time(tstop)
+         print*,'TOTAL elapsed time in thermals : ',tstop-tstart
 
 
@@ -156,49 +165,51 @@
 
 
-         DO l=2,nlayer
+         DO l=2,nlayermx
               pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))/ptimestep
          ENDDO
 
          pbl_dtke(:,1)=0.5*dq2_therm(:,1)/ptimestep
-         pbl_dtke(:,nlayer+1)=0.
+         pbl_dtke(:,nlayermx+1)=0.
 !! DIAGNOSTICS
         
         if(outptherm) then
-        if (ngrid .eq. 1) then
-        call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',&
+        if (ngridmx .eq. 1) then
+        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
      &                       'kg/m-2',1,entr_therm)
-        call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',&
+        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
      &                       'kg/m-2',1,detr_therm)
-        call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',&
+        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
      &                       'kg/m-2',1,fm_therm)
-        call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',&
+        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
      &                       'm/s',1,zw2)
-        call WRITEDIAGFI(ngrid,'heatFlux_up','heatFlux_updraft',&
+        call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
      &                       'SI',1,heatFlux)
-       call WRITEDIAGFI(ngrid,'heatFlux_down','heatFlux_downdraft',&
+       call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
      &                       'SI',1,heatFlux_down)
-        call WRITEDIAGFI(ngrid,'fraca','fraction coverage',&
+        call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
      &                       'percent',1,fraca)
-        call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',&
+        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
      &                       'm.s-2',1,buoyancyOut)
-        call WRITEDIAGFI(ngrid,'buoyancyEst','buoyancyEst',&
+        call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
      &                       'm.s-2',1,buoyancyEst)
-        call WRITEDIAGFI(ngrid,'d_t_th',  &
+        call WRITEDIAGFI(ngridmx,'d_t_th',  &
      &         'tendance temp TH','K',1,d_t_ajs)
+        call WRITEDIAGFI(ngridmx,'zmax',  &
+     &         'pbl height','m',0,zmax_th)
       else
 
-        call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',&
+        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
      &                       'kg/m-2',3,entr_therm)
-        call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',&
+        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
      &                       'kg/m-2',3,detr_therm)
-        call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',&
+        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
      &                       'kg/m-2',3,fm_therm)
-        call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',&
+        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
      &                       'm/s',3,zw2)
-        call WRITEDIAGFI(ngrid,'heatFlux','heatFlux',&
+        call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
      &                       'SI',3,heatFlux)
-        call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',&
+        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
      &                       'SI',3,buoyancyOut)
-        call WRITEDIAGFI(ngrid,'d_t_th',  &
+        call WRITEDIAGFI(ngridmx,'d_t_th',  &
      &         'tendance temp TH','K',3,d_t_ajs)
 
Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90	(revision 185)
@@ -2,50 +2,52 @@
 ! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $
 !
-      subroutine calltherm_mars(ngrid,nlayer,dtime,nq,zzlev,zzlay  &
+      subroutine calltherm_mars(dtime,zzlev,zzlay  &
      &      ,pplay,paprs,pphi  &
      &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
      &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,dq2_therm  &
-     &      ,fm_therm,entr_therm,detr_therm,lmax,&
+     &      ,fm_therm,entr_therm,detr_therm,lmax,zmax,&
      &   zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,&
      &     buoyancyOut,buoyancyEst,hfmax,wmax)
 
-       USE thermcell, only : nsplit_thermals,r_aspect_thermals
        USE ioipsl_getincom
       implicit none
 
-      INTEGER, INTENT(IN) :: ngrid,nlayer
+#include "dimensions.h"
+#include "dimphys.h"
+
       REAL dtime
-      LOGICAL logexpr0, logexpr2(ngrid,nlayer), logexpr1(ngrid)
+      LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx)
       REAL fact
-      INTEGER nbptspb,nq
-
-      REAL, INTENT(IN) :: zzlay(ngrid,nlayer)
-      REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1)
-
-      REAL u_seri(ngrid,nlayer),v_seri(ngrid,nlayer)
-      REAL t_seri(ngrid,nlayer),pq_therm(ngrid,nlayer,nq)
-      REAL q2_therm(ngrid,nlayer)
-      REAL paprs(ngrid,nlayer+1)
-      REAL pplay(ngrid,nlayer)
-      REAL pphi(ngrid,nlayer)
-      real zlev(ngrid,nlayer+1) 
+      INTEGER nbptspb
+
+      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
+
+      REAL u_seri(ngridmx,nlayermx),v_seri(ngridmx,nlayermx)
+      REAL t_seri(ngridmx,nlayermx),pq_therm(ngridmx,nlayermx,nqmx)
+      REAL q2_therm(ngridmx,nlayermx)
+      REAL paprs(ngridmx,nlayermx+1)
+      REAL pplay(ngridmx,nlayermx)
+      REAL pphi(ngridmx,nlayermx)
+      real zlev(ngridmx,nlayermx+1) 
 !test: on sort lentr et a* pour alimenter KE
-      REAL zw2(ngrid,nlayer+1),fraca(ngrid,nlayer+1)
-      REAL zzw2(ngrid,nlayer+1)
+      REAL zw2(ngridmx,nlayermx+1),fraca(ngridmx,nlayermx+1)
+      REAL zzw2(ngridmx,nlayermx+1)
 
 !FH Update Thermiques
-      REAL d_t_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
-      REAL d_u_ajs(ngrid,nlayer),d_v_ajs(ngrid,nlayer)
-      REAL dq2_therm(ngrid,nlayer), dq2_the(ngrid,nlayer)
-      real fm_therm(ngrid,nlayer+1)
-      real entr_therm(ngrid,nlayer),detr_therm(ngrid,nlayer)
+      REAL d_t_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
+      REAL d_u_ajs(ngridmx,nlayermx),d_v_ajs(ngridmx,nlayermx)
+      REAL dq2_therm(ngridmx,nlayermx), dq2_the(ngridmx,nlayermx)
+      real fm_therm(ngridmx,nlayermx+1)
+      real entr_therm(ngridmx,nlayermx),detr_therm(ngridmx,nlayermx)
 
 !********************************************************
 !     declarations
-      real zpopsk(ngrid,nlayer)
-      real ztla(ngrid,nlayer)
-      real wmax(ngrid)
-      real hfmax(ngrid)
-      integer lmax(ngrid)
+      real zpopsk(ngridmx,nlayermx)
+      real ztla(ngridmx,nlayermx)
+      real wmax(ngridmx)
+      real hfmax(ngridmx)
+      integer lmax(ngridmx)
+      real zmax(ngridmx)
 
 !nouvelles variables pour la convection
@@ -56,18 +58,20 @@
 
 ! variables locales
-      REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq)
-      REAL d_u_the(ngrid,nlayer),d_v_the(ngrid,nlayer)
+      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
+      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
 !
-      integer isplit
-      real zfm_therm(ngrid,nlayer+1),zdt
-      real zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer)
-      real heatFlux(ngrid,nlayer)
-      real heatFlux_down(ngrid,nlayer)
-      real buoyancyOut(ngrid,nlayer)
-      real buoyancyEst(ngrid,nlayer)
-      real zheatFlux(ngrid,nlayer)
-      real zheatFlux_down(ngrid,nlayer)
-      real zbuoyancyOut(ngrid,nlayer)
-      real zbuoyancyEst(ngrid,nlayer)
+      integer isplit,nsplit_thermals
+      real r_aspect_thermals
+
+      real zfm_therm(ngridmx,nlayermx+1),zdt
+      real zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
+      real heatFlux(ngridmx,nlayermx)
+      real heatFlux_down(ngridmx,nlayermx)
+      real buoyancyOut(ngridmx,nlayermx)
+      real buoyancyEst(ngridmx,nlayermx)
+      real zheatFlux(ngridmx,nlayermx)
+      real zheatFlux_down(ngridmx,nlayermx)
+      real zbuoyancyOut(ngridmx,nlayermx)
+      real zbuoyancyEst(ngridmx,nlayermx)
 
       character (len=20) :: modname='calltherm'
@@ -77,8 +81,12 @@
       logical, save :: first=.true.
 
+      REAL tstart,tstop
+
+
 !  Modele du thermique
 !  ===================
-
-         nsplit_thermals=20
+ 
+         r_aspect_thermals=3.
+         nsplit_thermals=40
          call getin("nsplit_thermals",nsplit_thermals)
 
@@ -98,4 +106,7 @@
          do isplit=1,nsplit_thermals
 
+!         call cpu_time(tstart)
+
+
 ! On reinitialise les flux de masse a zero pour le cumul en
 ! cas de splitting
@@ -116,13 +127,13 @@
          d_v_the(:,:)=0.
          dq2_the(:,:)=0.
-         if (nq .ne. 0) then
+         if (nqmx .ne. 0) then
             d_q_the(:,:,:)=0.
          endif
 
-             CALL thermcell_main_mars(ngrid,nlayer,nq,zdt  &
+             CALL thermcell_main_mars(zdt  &
      &      ,pplay,paprs,pphi,zzlev,zzlay  &
      &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
      &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
-     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax  &
+     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
      &      ,r_aspect_thermals &
      &      ,zzw2,fraca,zpopsk &
@@ -138,5 +149,5 @@
             dq2_the(:,:)=dq2_the(:,:)*fact            
 
-            if (nq .ne. 0) then
+            if (nqmx .ne. 0) then
                d_q_the(:,:,:)=d_q_the(:,:,:)*fact
             endif
@@ -175,4 +186,8 @@
             q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
 
+
+!         call cpu_time(tstop)
+!         print*,'elapsed time in thermals : ',tstop-tstart
+
          enddo ! isplit
 
@@ -180,6 +195,6 @@
 !****************************************************************
 
-!          do i=1,ngrid
-!             do k=1,nlayer
+!          do i=1,ngridmx
+!             do k=1,nlayermx
 !                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
 !               print*,'youpi je sers a quelque chose !'
@@ -187,5 +202,5 @@
 !          enddo
         
-          DO i=1,ngrid
+          DO i=1,ngridmx
             hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:))
             wmax(i)=MAXVAL(zw2(i,:))
Index: /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F	(revision 185)
@@ -361,10 +361,13 @@
 c Variables from thermal
 
-      REAL lmax_th_out(ngrid)
-      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
-      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
-      INTEGER lmax_th(ngrid)
-      REAL dtke_th(ngrid,nlayer+1)
-      REAL wmax_th(ngrid),hfmax_th(ngrid)
+      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
+      REAL wmax_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)
+      INTEGER lmax_th(ngridmx)
+      REAL dtke_th(ngridmx,nlayermx+1)
+      REAL dummycol(ngridmx)
+
 c=======================================================================
 #ifdef MESOSCALE
@@ -403,5 +406,5 @@
       PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngridmx,nsoilmx)
       mlayer(0:nsoilmx-1)=wdsoil(1,:)
-      PRINT*,'check: layer ', mlayer
+      PRINT*,'check: midlayer ', mlayer(:)
             !!!!!!!!!!!!!!!!! DONE in soil_setting.F 
             ! 1.5 Build layer(); following the same law as mlayer()
@@ -413,4 +416,7 @@
               layer(iloop)=lay1*(alpha**(iloop-1))
             enddo
+
+      PRINT*,'check: layer ', layer(:)
+
             !!!!!!!!!!!!!!!!! DONE in soil_setting.F
       tnom(:)=wtnom(:)   !! est rempli dans advtrac.h
@@ -1009,10 +1015,11 @@
       if(calltherm) then
  
-        call calltherm_interface(ngrid,nlayer,firstcall,
+        call calltherm_interface(firstcall,
      $ long,lati,zzlev,zzlay,
      $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
-     $ pplay,pplev,pphi,nq,zpopsk,
-     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th,hfmax_th,wmax_th)
- 
+     $ pplay,pplev,pphi,zpopsk,
+     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
+     $ dtke_th,hfmax_th,wmax_th) 
+
          DO l=1,nlayer
            DO ig=1,ngrid
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 185)
@@ -206,5 +206,6 @@
 
       CHARACTER*80 fichier 
-      INTEGER l,ig,ierr,igout,iq,i, tapphys
+      INTEGER l,ig,ierr,igout,iq,i, tapphys,k
+      LOGICAL end_of_file
 
       REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
@@ -302,12 +303,12 @@
 c Variables from thermal
 
-      REAL lmax_th_out(ngrid)
-      REAL wmax_th(ngrid),hfmax_th(ngrid)
-      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
-      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
-      INTEGER lmax_th(ngrid)
-      REAL dtke_th(ngrid,nlayer+1)
-      REAL dummycol(ngrid)
-      REAL hfx(ngrid)
+      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
+      REAL wmax_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)
+      INTEGER lmax_th(ngridmx)
+      REAL dtke_th(ngridmx,nlayermx+1)
+      REAL dummycol(ngridmx)
 
 c=======================================================================
@@ -562,6 +563,4 @@
            ENDIF
 
-
-
         ENDIF ! of if(mod(icount-1,iradia).eq.0)
 
@@ -578,5 +577,4 @@
                fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
            ENDDO
-
 
          DO l=1,nlayer
@@ -624,5 +622,5 @@
             enddo
          enddo
-         
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -634,6 +632,4 @@
      &        zdqdif,zdqsdif)
 
-          hfx(:) = zflubid(:)-capcal(:)*zdtsdif(:)
-
          DO l=1,nlayer
             DO ig=1,ngrid
@@ -647,7 +643,7 @@
          ENDDO
 
-         DO ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-         ENDDO
+          DO ig=1,ngrid
+             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+          ENDDO
 
          if (tracer) then 
@@ -680,9 +676,10 @@
       if(calltherm) then
  
-        call calltherm_interface(ngrid,nlayer,firstcall,
+        call calltherm_interface(firstcall,
      $ long,lati,zzlev,zzlay,
      $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
-     $ pplay,pplev,pphi,nq,zpopsk,
-     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th,hfmax_th,wmax_th)
+     $ pplay,pplev,pphi,zpopsk,
+     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
+     $ dtke_th,hfmax_th,wmax_th)
  
          DO l=1,nlayer
@@ -767,5 +764,5 @@
      $              co2ice,albedo,emis,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
-     $	            fluxsurf_sw,zls) 
+     $              fluxsurf_sw,zls) 
 
          DO l=1,nlayer
@@ -1482,8 +1479,4 @@
         lmax_th_out(:)=real(lmax_th(:))
 
-        call WRITEDIAGFI(ngridmx,'hfx',
-     &              'sensible heat flux','W.m-2',
-     &                         2,hfx)
-
         call WRITEDIAGFI(ngridmx,'lmax_th',
      &              'hauteur du thermique','K',
@@ -1534,7 +1527,4 @@
 ! THERMALS STUFF 1D
 
-        call WRITEDIAGFI(ngridmx,'hfx',
-     &              'sensible heat flux','W.m-2',
-     &                         0,hfx)
          if(calltherm) then
 
@@ -1580,8 +1570,4 @@
      &                  tsurf)
 
-           call WRITEDIAGFI(ngrid,"fluxsurf_sw",
-     &                "Solar radiative flux to surface","W.m-2",0,
-     &                fluxsurf_sw_tot)
-
          call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
          call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
Index: unk/LMDZ.MARS/libf/phymars/thermcell.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell.F90	(revision 184)
+++ 	(revision )
@@ -1,108 +1,0 @@
-MODULE thermcell
-!      USE ioipsl_getincom
-      IMPLICIT NONE
-
-      INTEGER :: iflag_thermals
-      REAL :: r_aspect_thermals,l_mix_thermals,tau_thermals
-      INTEGER :: w2di_thermals
-      INTEGER :: nsplit_thermals
-      INTEGER :: isplit
-      INTEGER :: iflag_coupl,iflag_clos,iflag_wake
-      INTEGER :: iflag_thermals_ed,iflag_thermals_optflux
-      REAL :: RG,RD
-      REAL :: rlvtt,rcpd,rtt,r2es
-      REAL :: retv,r5les,r5ies,rkappa
-      REAL :: R4LES,R4IES,R3IES,R3LES
-      INTEGER :: klon,klev
-      INTEGER :: prt_level,lunout
-      real,allocatable :: rlatd(:)
-      real,allocatable :: rlond(:)
-
-! Thermodynamic constants. The [OK] means that the constant has been made
-! compatible with the Martian model
-! ----------------------------------
-
-! RG : mars mean gravity field                          [OK]
-! RD : dry air constant = 1000 R/Mair                   [OK]
-! rlvtt : value of Lv(Tt) (vapo. latent heat at Tt)     [dep. on air ?]
-! rcpd : Cp of dry air = 7/2 RD for a perfect gas       [OK,CHECK VALUE with gcm one]
-!rcpd = r/(r/CP), avec r/CP = 1/4.4 pour LES, 0.256793 pour gcm (1./3.89419), et r=1000.R/Mair=RD
-! rtt : triple point temperature Tt                     [OK]
-! r2es :                                                [Probably OK]
-! retv = restt*RD/RV :                                  [-]
-! restt : saturation pressure at Tt es(Tt) = 611.14 Pa  [dep. on air ?]
-! RV : vapor constant = 1000 R/Mh2o                     [OK]
-! r3les :                                               [Probably OK]
-! r3ies :                                               [Probably OK]
-! r4les :                                               [Probably OK]
-! r4ies :                                               [Probably OK]
-! r5ies = r3ies*(rtt-r4les)                             [Probably OK]
-! rkappa = RD/RCPD = r/cp = 0.256793  			[OK]
-! FOEEW : vapeur d'eau saturante			[OK]
-! FOEDE : derive par rapport a la temperature		[OK]
-      PARAMETER (RG=3.72,RD=191.182)
-      PARAMETER (rlvtt=2.5008E+06,rcpd=RD*3.89419)
-      PARAMETER (rtt=273.16,r2es=253.156)
-      PARAMETER (retv=1.41409,r5les=4097.93)
-      PARAMETER (r5ies=5807.81,rkappa=1./3.89419)
-      PARAMETER (R4IES=7.66,R4LES=35.86)
-      PARAMETER (R3IES=21.875,R3LES=17.269)
-
-! r_aspect_thermals : rapport d'aspect : parameter defined in FH's HDR
-!   length over height ratio of the thermals in the alimentation layer (check)
-!   it's value is advised to be set to 2 in HDR           [OK]
-! l_mix_thermals : mixing length => lambda : parameter for plume penetration in
-! the atm (the plum size will decrease with sqrt(lambda) : USELESS
-! w2di_thermals : ?
-! tau_thermals : relaxation length of the thermals : the model implements the thermals 
-!   tendancies with a typical relaxation time. Typical value is said to be 1800 s in the
-!   code however, it is found at 720 s in the .def files  [OK]
-! nsplit_thermals : subdivisions of the timestep in the calculations of the plume.
-!   Increase the calculation precision greatly, but also computation time
-!   typical value is found at 1 in .def (no split)        [OK]
-! prt_level : print level for the gcm printed outputs     [OK]
-! lunout : output channel for prt level                   [OK]
-! iflag_thermals : choice of the thermals version. 15 and 16 are the newest ones
-!   and 16 activates "bidouille stratocu" which is required for now
-! 							  [OK]
-! iflag_thermals_ed : choice of the plume version
-!   8 is the normal version, 9 is the working version of AJ, 10 is CR
-!   also used for cases .ge.1 in thermcell_height         [OK]
-! iflag_coupl : USELESS ? 				  [useless ?]
-! iflag_clos : USELESS ?				  [useless ?]
-! iflag_wake : USELESS ?				  [useless ?]
-      PARAMETER (r_aspect_thermals = 1.5)
-      PARAMETER (l_mix_thermals = 30.)
-      PARAMETER (w2di_thermals = 1)
-!      PARAMETER (tau_thermals = 720.)
-!      PARAMETER (nsplit_thermals = 1)
-      PARAMETER (prt_level=0,lunout=6)
-      PARAMETER (iflag_thermals=20,iflag_thermals_ed=10)
-      PARAMETER (iflag_thermals_optflux=1,iflag_coupl=1)
-      PARAMETER (iflag_clos=2)
-
-      CONTAINS
-        FUNCTION foede(PTARG,PDELARG,P5ARG,PQSARG,PCOARG)
-          IMPLICIT NONE
-          REAL :: foede
-          REAL, INTENT(IN) :: PTARG,PDELARG,P5ARG,PQSARG,PCOARG
-          foede = PQSARG*PCOARG*P5ARG &
-       & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
-        END FUNCTION foede
-
-        FUNCTION foeew(PTARG,PDELARG)
-          IMPLICIT NONE
-          REAL :: foeew
-          REAL, INTENT(IN) :: PTARG,PDELARG
-          foeew = EXP (                                   &
-       &          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)        &
-       & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
-
-        END FUNCTION foeew
-
-! Variables which have been moved to .def and called independantly :
-! nsplit_thermals, tau
-
-END MODULE thermcell
-
-
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90	(revision 185)
@@ -1,3 +1,3 @@
-      subroutine thermcell_dqupdown(ngrid,nlay,ptimestep,fm0,entr0,  &
+      subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0,  &
      &    detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
       implicit none
@@ -13,34 +13,37 @@
 !=======================================================================
 
+#include "dimensions.h"
+#include "dimphys.h"
+
 ! ============================ INPUTS ============================
 
-      INTEGER, INTENT(IN) :: ngrid,nlay
+      INTEGER, INTENT(IN) :: ngrid,nlayer
       REAL, INTENT(IN) :: ptimestep
-      REAL, INTENT(IN) :: fm0(ngrid,nlay+1)
-      REAL, INTENT(IN) :: entr0(ngrid,nlay),detr0(ngrid,nlay)
-      REAL, INTENT(IN) :: q_therm(ngrid,nlay)
-      REAL, INTENT(IN) :: fm_down(ngrid,nlay+1)
-      REAL, INTENT(IN) :: ztvd(ngrid,nlay)
-      REAL, INTENT(IN) :: ztv(ngrid,nlay)
+      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
+      REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1)
+      REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: ztv(ngridmx,nlayermx)
       CHARACTER (LEN=20), INTENT(IN) :: charvar
-      REAL, INTENT(IN) :: masse0(ngrid,nlay)
-      INTEGER, INTENT(IN) :: lmax(ngrid)
+      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
+      INTEGER, INTENT(IN) :: lmax(ngridmx)
 
 ! ============================ OUTPUTS ===========================
 
-      REAL, INTENT(OUT) :: dq_therm(ngrid,nlay)
+      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)
 
 ! ============================ LOCAL =============================
 
-!      REAL detr0(ngrid,nlay)
-      REAL detrd(ngrid,nlay)
-      REAL entrd(ngrid,nlay)      
-      REAL fmd(ngrid,nlay+1)
-      REAL q(ngrid,nlay)
-      REAL qa(ngrid,nlay)
-      REAL qd(ngrid,nlay)
+!      REAL detr0(ngridmx,nlayermx)
+      REAL detrd(ngridmx,nlayermx)
+      REAL entrd(ngridmx,nlayermx)      
+      REAL fmd(ngridmx,nlayermx+1)
+      REAL q(ngridmx,nlayermx)
+      REAL qa(ngridmx,nlayermx)
+      REAL qd(ngridmx,nlayermx)
       INTEGER ig,k
-      LOGICAL active(ngrid,nlay)
-      INTEGER lmax_down(ngrid),lmin_down(ngrid)
+      LOGICAL active(ngridmx,nlayermx)
+      INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
       INTEGER ncorec
 
@@ -64,14 +67,14 @@
 !! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
 !
-!      do ig=1,ngrid
-!          if (ztv(ig,nlay)-ztvd(ig,nlay) .gt. 0.5) then
+!      do ig=1,ngridmx
+!          if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
 !            print*,"downdraft non nul derniere couche !!! (dqupdown)"
 !          endif
-!          detrd(ig,nlay)=0.
-!          entrd(ig,nlay)=0.
-!      enddo
-!
-!      do k=nlay-1,1,-1
-!          do ig=1,ngrid
+!          detrd(ig,nlayermx)=0.
+!          entrd(ig,nlayermx)=0.
+!      enddo
+!
+!      do k=nlayermx-1,1,-1
+!          do ig=1,ngridmx
 !
 !          if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
@@ -93,6 +96,6 @@
 !      lmin_down(:)=1
 !
-!      do k=1,nlay
-!        do ig=1,ngrid
+!      do k=1,nlayermx
+!        do ig=1,ngridmx
 !        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
 !!         if (entrd(ig,k).gt.detrd(ig,k)) then
@@ -101,6 +104,6 @@
 !        enddo
 !      enddo
-!      do k=nlay,1,-1
-!        do ig=1,ngrid
+!      do k=nlayermx,1,-1
+!        do ig=1,ngridmx
 !        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
 !!         if (detrd(ig,k).gt.entrd(ig,k)) then
@@ -113,5 +116,5 @@
 !      fmd(:,:)=0.
 !       
-!      do ig=1,ngrid
+!      do ig=1,ngridmx
 !          if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
 !!          fmd(ig,lmax_down(ig))=0.
@@ -137,6 +140,6 @@
 !      enddo
 !         ncorec=0
-!         do k=nlay,2,-1
-!         do ig=1,ngrid
+!         do k=nlayermx,2,-1
+!         do ig=1,ngridmx
 !            if (fmd(ig,k).lt.0.) then
 !!               detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
@@ -160,6 +163,6 @@
 !         print*, lmin_down(:),lmax_down(:)
 !
-!      do k=2,nlay
-!      do ig=1,ngrid
+!      do k=2,nlayermx
+!      do ig=1,ngridmx
 !          active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
 !     & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k))
@@ -168,5 +171,5 @@
 !
 !
-!      do ig=1,ngrid
+!      do ig=1,ngridmx
 !      do k=lmin_down(ig),lmax_down(ig)
 !          if(.not.active(ig,k)) then
@@ -180,10 +183,10 @@
 !      endif
 !
-!!      do ig=1,ngrid
+!!      do ig=1,ngridmx
 !!         active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
 !!     &       ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
 !!      enddo
 !!
-!!      do ig=1,ngrid
+!!      do ig=1,ngridmx
 !!          if (lmax_down(ig).gt.1) then
 !!            do k=lmax_down(ig)-1,lmin_down(ig),-1
@@ -199,6 +202,6 @@
 !! ========== qa : q in updraft ==================================
 !
-      do k=2,nlay
-         do ig=1,ngrid
+      do k=2,nlayermx
+         do ig=1,ngridmx
             if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt.  &
      &         1.e-5*masse0(ig,k)) then
@@ -219,6 +222,6 @@
 ! ========== qd : q in downdraft =================================
 !
-!      do k=nlay-1,1,-1
-!         do ig=1,ngrid
+!      do k=nlayermx-1,1,-1
+!         do ig=1,ngridmx
 !            if (active(ig,k)) then
 !         qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k))  &
@@ -241,5 +244,5 @@
 ! ====== dq ======================================================
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if(active(ig,1)) then
 
@@ -258,6 +261,6 @@
        enddo
       
-       do k=2,nlay-1
-         do ig=1, ngrid
+       do k=2,nlayermx-1
+         do ig=1, ngridmx
 
          if(active(ig,k)) then
@@ -281,18 +284,18 @@
       enddo
 
-         do ig=1, ngrid
-
-         if(active(ig,nlay)) then
-
-         dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay)+detrd(ig,nlay)*qd(ig,nlay) &
-      &               +fmd(ig,nlay)*q(ig,nlay-1)   &
-      &         -entr0(ig,nlay)*q(ig,nlay)-entrd(ig,nlay)*q(ig,nlay)   &
-      &               -fm0(ig,nlay)*q(ig,nlay)) &
-      &               *ptimestep/masse0(ig,nlay)
+         do ig=1, ngridmx
+
+         if(active(ig,nlayermx)) then
+
+         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx)+detrd(ig,nlayermx)*qd(ig,nlayermx) &
+      &               +fmd(ig,nlayermx)*q(ig,nlayermx-1)   &
+      &         -entr0(ig,nlayermx)*q(ig,nlayermx)-entrd(ig,nlayermx)*q(ig,nlayermx)   &
+      &               -fm0(ig,nlayermx)*q(ig,nlayermx)) &
+      &               *ptimestep/masse0(ig,nlayermx)
 
          else
-         dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay) &
-      &             -entr0(ig,nlay)*q(ig,nlay)-fm0(ig,nlay)*q(ig,nlay)) &
-      &               *ptimestep/masse0(ig,nlay)
+         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) &
+      &             -entr0(ig,nlayermx)*q(ig,nlayermx)-fm0(ig,nlayermx)*q(ig,nlayermx)) &
+      &               *ptimestep/masse0(ig,nlayermx)
          endif
          
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 184)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 185)
@@ -1,9 +1,9 @@
 !
 !
-      SUBROUTINE thermcell_main_mars(ngrid,nlay,nq,ptimestep  &
+      SUBROUTINE thermcell_main_mars(ptimestep  &
      &                  ,pplay,pplev,pphi,zlev,zlay  &
      &                  ,pu,pv,pt,pq,pq2  &
      &                  ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj  &
-     &                  ,fm,entr,detr,lmax  &
+     &                  ,fm,entr,detr,lmax,zmax  &
      &                  ,r_aspect &
      &                  ,zw2,fraca &
@@ -11,5 +11,4 @@
      &                  ,buoyancyOut, buoyancyEst)
 
-      USE thermcell,ONLY:RG,RD,RKAPPA
       IMPLICIT NONE
 
@@ -17,40 +16,44 @@
 ! Mars version of thermcell_main. Author : A Colaitis
 !=======================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+
 ! ============== INPUTS ==============
 
-      INTEGER, INTENT(IN) :: ngrid,nlay,nq
       REAL, INTENT(IN) :: ptimestep,r_aspect
-      REAL, INTENT(IN) :: pt(ngrid,nlay)
-      REAL, INTENT(IN) :: pu(ngrid,nlay)
-      REAL, INTENT(IN) :: pv(ngrid,nlay)
-      REAL, INTENT(IN) :: pq(ngrid,nlay,nq)
-      REAL, INTENT(IN) :: pq2(ngrid,nlay)
-      REAL, INTENT(IN) :: pplay(ngrid,nlay)
-      REAL, INTENT(IN) :: pplev(ngrid,nlay+1)
-      REAL, INTENT(IN) :: pphi(ngrid,nlay)
-      REAL, INTENT(IN) :: zlay(ngrid,nlay)
-      REAL, INTENT(IN) :: zlev(ngrid,nlay+1)
+      REAL, INTENT(IN) :: pt(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pu(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pv(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx)
+      REAL, INTENT(IN) :: pq2(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
+      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: zlay(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1)
 
 ! ============== OUTPUTS ==============
 
-      REAL, INTENT(OUT) :: pdtadj(ngrid,nlay)
-      REAL, INTENT(OUT) :: pduadj(ngrid,nlay)
-      REAL, INTENT(OUT) :: pdvadj(ngrid,nlay)
-      REAL, INTENT(OUT) :: pdqadj(ngrid,nlay,nq)
-      REAL, INTENT(OUT) :: pdq2adj(ngrid,nlay)
-      REAL, INTENT(OUT) :: zw2(ngrid,nlay+1)
+      REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
+      REAL, INTENT(OUT) :: pduadj(ngridmx,nlayermx)
+      REAL, INTENT(OUT) :: pdvadj(ngridmx,nlayermx)
+      REAL, INTENT(OUT) :: pdqadj(ngridmx,nlayermx,nqmx)
+      REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
+      REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
 
 ! Diagnostics
-      REAL, INTENT(OUT) :: heatFlux(ngrid,nlay)   ! interface heatflux
-      REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft
-      REAL, INTENT(OUT) :: buoyancyOut(ngrid,nlay)  ! interlayer buoyancy term
-      REAL, INTENT(OUT) :: buoyancyEst(ngrid,nlay)  ! interlayer estimated buoyancy term
+      REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
+      REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
+      REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
+      REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
 
 ! dummy variables when output not needed :
 
-!      REAL :: heatFlux(ngrid,nlay)   ! interface heatflux
-!      REAL :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft
-!      REAL :: buoyancyOut(ngrid,nlay)  ! interlayer buoyancy term
-!      REAL :: buoyancyEst(ngrid,nlay)  ! interlayer estimated buoyancy term
+!      REAL :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
+!      REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
+!      REAL :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
+!      REAL :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
 
 
@@ -58,45 +61,45 @@
 
       INTEGER ig,k,l,ll,iq
-      INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
-      INTEGER lmix(ngrid)
-      INTEGER lmix_bis(ngrid)
-      REAL linter(ngrid)
-      REAL zmix(ngrid)
-      REAL zmax(ngrid)
-      REAL ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay)
-      REAL zmax_sec(ngrid)
-      REAL zh(ngrid,nlay)
-      REAL zdthladj(ngrid,nlay)
-      REAL zdthladj_down(ngrid,nlay)
-      REAL ztvd(ngrid,nlay)
-      REAL ztv(ngrid,nlay)
-      REAL zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay)
-      REAL zva(ngrid,nlay)
-      REAL zua(ngrid,nlay)
-
-      REAL zta(ngrid,nlay)
-      REAL fraca(ngrid,nlay+1)
-      REAL q2(ngrid,nlay)
-      REAL rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay)
-      REAL zpopsk(ngrid,nlay)
-
-      REAL wmax(ngrid)
-      REAL wmax_sec(ngrid)
-      REAL fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay)
-
-      REAL fm_down(ngrid,nlay+1)
-
-      REAL ztla(ngrid,nlay)
-
-      REAL f_star(ngrid,nlay+1),entr_star(ngrid,nlay)
-      REAL detr_star(ngrid,nlay)
-      REAL alim_star_tot(ngrid)
-      REAL alim_star(ngrid,nlay)
-      REAL alim_star_clos(ngrid,nlay)
-      REAL f(ngrid)
-
-      REAL teta_th_int(ngrid,nlay)
-      REAL teta_env_int(ngrid,nlay)
-      REAL teta_down_int(ngrid,nlay)
+      INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
+      INTEGER lmix(ngridmx)
+      INTEGER lmix_bis(ngridmx)
+      REAL linter(ngridmx)
+      REAL zmix(ngridmx)
+      REAL zmax(ngridmx)
+      REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
+      REAL zmax_sec(ngridmx)
+      REAL zh(ngridmx,nlayermx)
+      REAL zdthladj(ngridmx,nlayermx)
+      REAL zdthladj_down(ngridmx,nlayermx)
+      REAL ztvd(ngridmx,nlayermx)
+      REAL ztv(ngridmx,nlayermx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx),zo(ngridmx,nlayermx)
+      REAL zva(ngridmx,nlayermx)
+      REAL zua(ngridmx,nlayermx)
+
+      REAL zta(ngridmx,nlayermx)
+      REAL fraca(ngridmx,nlayermx+1)
+      REAL q2(ngridmx,nlayermx)
+      REAL rho(ngridmx,nlayermx),rhobarz(ngridmx,nlayermx),masse(ngridmx,nlayermx)
+      REAL zpopsk(ngridmx,nlayermx)
+
+      REAL wmax(ngridmx)
+      REAL wmax_sec(ngridmx)
+      REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
+
+      REAL fm_down(ngridmx,nlayermx+1)
+
+      REAL ztla(ngridmx,nlayermx)
+
+      REAL f_star(ngridmx,nlayermx+1),entr_star(ngridmx,nlayermx)
+      REAL detr_star(ngridmx,nlayermx)
+      REAL alim_star_tot(ngridmx)
+      REAL alim_star(ngridmx,nlayermx)
+      REAL alim_star_clos(ngridmx,nlayermx)
+      REAL f(ngridmx)
+
+      REAL teta_th_int(ngridmx,nlayermx)
+      REAL teta_env_int(ngridmx,nlayermx)
+      REAL teta_down_int(ngridmx,nlayermx)
 
       CHARACTER (LEN=20) :: modname
@@ -105,9 +108,10 @@
 ! ============= PLUME VARIABLES ============
 
-      REAL w_est(ngrid,nlay+1)
-      REAL wa_moy(ngrid,nlay+1)
-      REAL wmaxa(ngrid)
-      REAL zdz,zbuoy(ngrid,nlay),zw2m
-      LOGICAL active(ngrid),activetmp(ngrid)
+      REAL w_est(ngridmx,nlayermx+1)
+      REAL wa_moy(ngridmx,nlayermx+1)
+      REAL wmaxa(ngridmx)
+      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
+      LOGICAL active(ngridmx),activetmp(ngridmx)
+      REAL a1,b1,ae,be,ad,bd
 
 ! ==========================================
@@ -115,7 +119,7 @@
 ! ============= HEIGHT VARIABLES ===========
 
-      REAL num(ngrid)
-      REAL denom(ngrid)
-      REAL zlevinter(ngrid)
+      REAL num(ngridmx)
+      REAL denom(ngridmx)
+      REAL zlevinter(ngridmx)
 
 ! =========================================
@@ -123,11 +127,11 @@
 ! ============= DRY VARIABLES =============
 
-       REAL zw2_dry(ngrid,nlay+1)
-       REAL f_star_dry(ngrid,nlay+1)
-       REAL ztva_dry(ngrid,nlay+1)
-       REAL wmaxa_dry(ngrid)
-       REAL wa_moy_dry(ngrid,nlay+1)
-       REAL linter_dry(ngrid),zlevinter_dry(ngrid)
-       INTEGER lmix_dry(ngrid),lmax_dry(ngrid)
+       REAL zw2_dry(ngridmx,nlayermx+1)
+       REAL f_star_dry(ngridmx,nlayermx+1)
+       REAL ztva_dry(ngridmx,nlayermx+1)
+       REAL wmaxa_dry(ngridmx)
+       REAL wa_moy_dry(ngridmx,nlayermx+1)
+       REAL linter_dry(ngridmx),zlevinter_dry(ngridmx)
+       INTEGER lmix_dry(ngridmx),lmax_dry(ngridmx)
 
 ! =========================================
@@ -135,7 +139,7 @@
 ! ============= CLOSURE VARIABLES =========
 
-      REAL zdenom(ngrid)
-      REAL alim_star2(ngrid)
-      REAL alim_star_tot_clos(ngrid)
+      REAL zdenom(ngridmx)
+      REAL alim_star2(ngridmx)
+      REAL alim_star_tot_clos(ngridmx)
       INTEGER llmax
 
@@ -161,10 +165,10 @@
 !               not used for now
 
-      REAL qa(ngrid,nlay),detr_dv2(ngrid,nlay),zf,zf2
-      REAL wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
-      REAL gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
-      REAL ue(ngrid,nlay),ve(ngrid,nlay)
-      LOGICAL ltherm(ngrid,nlay)
-      REAL dua(ngrid,nlay),dva(ngrid,nlay)
+      REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2
+      REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)
+      REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)
+      REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)
+      LOGICAL ltherm(ngridmx,nlayermx)
+      REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)
       INTEGER iter
       INTEGER nlarga0
@@ -202,26 +206,26 @@
 !-----------------------------------------------------------------------
 
-!      do l=2,nlay
-!         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
+!      do l=2,nlayermx
+!         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
 !      enddo
 !         zlev(:,1)=0.
-!         zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
-
-!         zlay(:,:)=pphi(:,:)/RG
+!         zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
+
+!         zlay(:,:)=pphi(:,:)/g
 !-----------------------------------------------------------------------
 !   Calcul des densites
 !-----------------------------------------------------------------------
 
-      rho(:,:)=pplay(:,:)/(RD*pt(:,:))
+      rho(:,:)=pplay(:,:)/(r*pt(:,:))
 
       rhobarz(:,1)=rho(:,1)
 
-      do l=2,nlay
+      do l=2,nlayermx
          rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
       enddo
 
 !calcul de la masse
-      do l=1,nlay
-         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
+      do l=1,nlayermx
+         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
       enddo
 
@@ -309,4 +313,10 @@
       wa_moy(:,:)=0.
       linter(:)=1.
+
+!      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
+
 ! Initialisation des variables entieres
       lmix(:)=1
@@ -320,5 +330,5 @@
 !-------------------------------------------------------------------------
       active(:)=ztv(:,1)>ztv(:,2)
-          do ig=1,ngrid
+          do ig=1,ngridmx
             if (ztv(ig,1)>=(ztv(ig,2))) then
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
@@ -330,7 +340,8 @@
          enddo
 
-       do l=2,nlay-1
-         do ig=1,ngrid
-            if (ztv(ig,l)>(ztv(ig,l+1)+0.5) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
+       do l=2,nlayermx-1
+!        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
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
      &                       *sqrt(zlev(ig,l+1))
@@ -341,6 +352,6 @@
         enddo
       enddo
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
             if (alim_star_tot(ig) > 1.e-10 ) then
                alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
@@ -350,7 +361,7 @@
 
       alim_star_tot(:)=1.
-      if(alim_star(1,1) .ne. 0.) then
-      print*, alim_star(:,:)
-      endif
+!      if(alim_star(1,1) .ne. 0.) then
+!      print*, 'lalim star' ,lalim(1)
+!      endif
 
 !------------------------------------------------------------------------------
@@ -362,5 +373,5 @@
 !------------------------------------------------------------------------------
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
 ! Le panache va prendre au debut les caracteristiques de l'air contenu
 ! dans cette couche.
@@ -371,5 +382,5 @@
           f_star(ig,1)=0.
           f_star(ig,2)=alim_star(ig,1)
-          zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
+          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))
@@ -383,5 +394,5 @@
 !boucle de calcul de la vitesse verticale dans le thermique
 !==============================================================================
-      do l=2,nlay-1
+      do l=2,nlayermx-1
 !==============================================================================
 
@@ -389,5 +400,5 @@
 ! On decide si le thermique est encore actif ou non
 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
-          do ig=1,ngrid
+          do ig=1,ngridmx
              active(ig)=active(ig) &
       &                 .and. zw2(ig,l)>1.e-10 &
@@ -405,5 +416,5 @@
 !---------------------------------------------------------------------------
 
-          do ig=1,ngrid
+          do ig=1,ngridmx
              if(active(ig)) then
 
@@ -417,10 +428,11 @@
 
                 zdz=zlev(ig,l+1)-zlev(ig,l)
-                zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-                if (((2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
-                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015 &
-     & -2.*zdz*w_est(ig,l)*0.045*(2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015)**0.6)
+                zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+
+                if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
+                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 &
+     & -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*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015)
+                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)
                 endif
                 if (w_est(ig,l+1).lt.0.) then
@@ -434,21 +446,21 @@
 !-------------------------------------------------
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
         if (active(ig)) then
 
           zw2m=w_est(ig,l+1)
-          if((2.5*(zbuoy(ig,l)/zw2m)-0.0015) .gt. 0.) then
+
+          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
-!       &   0.0118*((zbuoy(ig,l)/zw2m)/0.043)**(1./1.65)
-!       &   0.0090*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.90)
-!       &   0.0120*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.6)
-        &   MAX(0.,0.045*(2.5*(zbuoy(ig,l)/zw2m)-0.0015)**0.6)
+        &   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.
+                detr_star(ig,l)=0.
              else
+
 !                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
 !              &     0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7)
@@ -462,5 +474,8 @@
 ! new param from continuity eq with a fit on dfdz
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
-               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)
+            &  ad
+
+!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
+!                &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)     
 
 !              &     0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35)
@@ -471,5 +486,7 @@
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
-               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)
+                &     bd*zbuoy(ig,l)/zw2m
+
+!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
 
 !              &  *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m)
@@ -506,5 +523,5 @@
 
       activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
-      do ig=1,ngrid
+      do ig=1,ngridmx
        if (activetmp(ig)) then
 
@@ -516,14 +533,14 @@
       enddo
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
       if (activetmp(ig)) then
            ztva(ig,l) = ztla(ig,l)
-           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-
-           if (((2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
-           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-         &
-     & 2.*zdz*zw2(ig,l)*0.0015-2.*zdz*zw2(ig,l)*0.045*(2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015)**0.6)
+           zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+
+           if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
+           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-         &
+     & 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*2.5*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*0.0015)
+           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1)
            endif
       endif
@@ -534,5 +551,5 @@
 !---------------------------------------------------------------------------
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
              print*,'On tombe sur le cas particulier de thermcell_plume'
@@ -562,8 +579,8 @@
 
 !on recalcule alim_star_tot
-       do ig=1,ngrid
+       do ig=1,ngridmx
           alim_star_tot(ig)=0.
        enddo
-       do ig=1,ngrid
+       do ig=1,ngridmx
           do l=1,lalim(ig)-1
           alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
@@ -571,6 +588,6 @@
        enddo
 
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
             if (alim_star_tot(ig) > 1.e-10 ) then
                alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
@@ -595,9 +612,9 @@
 
 !calcul de la hauteur max du thermique
-      do ig=1,ngrid
+      do ig=1,ngridmx
          lmax(ig)=lalim(ig)
       enddo
-      do ig=1,ngrid
-         do l=nlay,lalim(ig)+1,-1
+      do ig=1,ngridmx
+         do l=nlayermx,lalim(ig)+1,-1
             if (zw2(ig,l).le.1.e-10) then
                lmax(ig)=l-1
@@ -608,13 +625,13 @@
 ! On traite le cas particulier qu'il faudrait éviter ou le thermique
 ! atteind le haut du modele ...
-      do ig=1,ngrid
-      if ( zw2(ig,nlay) > 1.e-10 ) then
+      do ig=1,ngridmx
+      if ( zw2(ig,nlayermx) > 1.e-10 ) then
           print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
-          lmax(ig)=nlay
+          lmax(ig)=nlayermx
       endif
       enddo
 
 ! pas de thermique si couche 1 stable
-!      do ig=1,ngrid
+!      do ig=1,ngridmx
 !         if (lmin(ig).gt.1) then
 !             lmax(ig)=1
@@ -625,10 +642,10 @@
 !
 ! Determination de zw2 max
-      do ig=1,ngrid
+      do ig=1,ngridmx
          wmax(ig)=0.
       enddo
 
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
             if (l.le.lmax(ig)) then
                 if (zw2(ig,l).lt.0.)then
@@ -643,5 +660,5 @@
       enddo
 !   Longueur caracteristique correspondant a la hauteur des thermiques.
-      do  ig=1,ngrid
+      do  ig=1,ngridmx
          zmax(ig)=0.
          zlevinter(ig)=zlev(ig,1)
@@ -650,11 +667,11 @@
          num(:)=0.
          denom(:)=0.
-         do ig=1,ngrid
-          do l=1,nlay
+         do ig=1,ngridmx
+          do l=1,nlayermx
              num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
              denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
           enddo
        enddo
-       do ig=1,ngrid
+       do ig=1,ngridmx
        if (denom(ig).gt.1.e-10) then
           zmax(ig)=2.*num(ig)/denom(ig)
@@ -663,5 +680,5 @@
 
 ! def de  zmix continu (profil parabolique des vitesses)
-      do ig=1,ngrid
+      do ig=1,ngridmx
            if (lmix(ig).gt.1) then
 ! test
@@ -673,7 +690,7 @@
 !
             zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
-     &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)  &
+     &        *((zlev(ig,lmix(ig))*zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)*zlev(ig,lmix(ig)+1)))  &
      &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
-     &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))  &
+     &        *((zlev(ig,lmix(ig)-1)*zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))*zlev(ig,lmix(ig)))))  &
      &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
      &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
@@ -694,6 +711,6 @@
 !
 ! calcul du nouveau lmix correspondant
-      do ig=1,ngrid
-         do l=1,nlay
+      do ig=1,ngridmx
+         do l=1,nlayermx
             if (zmix(ig).ge.zlev(ig,l).and.  &
      &          zmix(ig).lt.zlev(ig,l+1)) then
@@ -710,132 +727,4 @@
 ! ===========================================================================
 
-
-! ===========================================================================
-! ================= DRY =====================================================
-! ===========================================================================
-
-!initialisations
-       do ig=1,ngrid
-          do l=1,nlay+1
-             zw2_dry(ig,l)=0.
-             wa_moy_dry(ig,l)=0.
-          enddo
-       enddo
-       do ig=1,ngrid
-          do l=1,nlay
-             ztva_dry(ig,l)=ztv(ig,l)
-          enddo
-       enddo
-       do ig=1,ngrid
-          wmax_sec(ig)=0.
-          wmaxa_dry(ig)=0.
-       enddo
-!calcul de la vitesse a partir de la CAPE en melangeant thetav
-
-
-! Calcul des F^*, integrale verticale de E^*
-       f_star_dry(:,1)=0.
-       do l=1,nlay
-          f_star_dry(:,l+1)=f_star_dry(:,l)+alim_star(:,l)
-       enddo
-
-! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
-       linter_dry(:)=0.
-
-! couche la plus haute concernee par le thermique.
-       lmax_dry(:)=1
-
-! Le niveau linter est une variable continue qui se trouve dans la couche
-! lmax
-
-       do l=1,nlay-2
-         do ig=1,ngrid
-            if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
-
-!------------------------------------------------------------------------
-!  Calcul de la vitesse en haut de la premiere couche instable.
-!  Premiere couche du panache thermique
-!------------------------------------------------------------------------
-
-             zw2_dry(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
-     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
-     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
-!------------------------------------------------------------------------
-! Tant que la vitesse en bas de la couche et la somme du flux de masse
-! et de l'entrainement (c'est a dire le flux de masse en haut) sont
-! positifs, on calcul
-! 1. le flux de masse en haut  f_star(ig,l+1)
-! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
-! 3. la vitesse au carré en haut zw2(ig,l+1)
-!------------------------------------------------------------------------
-
-            else if (zw2_dry(ig,l).ge.1e-10) then
-
-       ztva_dry(ig,l)=(f_star_dry(ig,l)*ztva_dry(ig,l-1)+alim_star(ig,l) &
-     &                    *ztv(ig,l))/f_star_dry(ig,l+1)
-      zw2_dry(ig,l+1)=zw2_dry(ig,l)*(f_star_dry(ig,l)/f_star_dry(ig,l+1))**2+  &
-     &                     2.*RG*(ztva_dry(ig,l)-ztv(ig,l))/ztv(ig,l)  &
-     &                     *(zlev(ig,l+1)-zlev(ig,l))
-            endif
-! determination de zmax continu par interpolation lineaire
-!------------------------------------------------------------------------
-
-            if (zw2_dry(ig,l+1)>0. .and. zw2_dry(ig,l+1).lt.1.e-10) then
-!               stop'On tombe sur le cas particulier de thermcell_dry'
-!               print*,'On tombe sur le cas particulier de thermcell_dry'
-                zw2_dry(ig,l+1)=0.
-                linter_dry(ig)=l+1
-                lmax_dry(ig)=l
-            endif
-
-            if (zw2_dry(ig,l+1).lt.0.) then
-               linter_dry(ig)=(l*(zw2_dry(ig,l+1)-zw2_dry(ig,l))  &
-     &           -zw2_dry(ig,l))/(zw2_dry(ig,l+1)-zw2_dry(ig,l))
-               zw2_dry(ig,l+1)=0.
-               lmax_dry(ig)=l
-            endif
-
-               wa_moy_dry(ig,l+1)=sqrt(zw2_dry(ig,l+1))
-
-            if (wa_moy_dry(ig,l+1).gt.wmaxa_dry(ig)) then
-!   lmix est le niveau de la couche ou w (wa_moy) est maximum
-               lmix_dry(ig)=l+1
-               wmaxa_dry(ig)=wa_moy_dry(ig,l+1)
-            endif
-         enddo
-      enddo
-!
-! Determination de zw2 max
-      do ig=1,ngrid
-         wmax_sec(ig)=0.
-      enddo
-      do l=1,nlay
-         do ig=1,ngrid
-            if (l.le.lmax_dry(ig)) then
-                zw2_dry(ig,l)=sqrt(zw2_dry(ig,l))
-                wmax_sec(ig)=max(wmax_sec(ig),zw2_dry(ig,l))
-            else
-                 zw2_dry(ig,l)=0.
-            endif
-          enddo
-      enddo
-
-!   Longueur caracteristique correspondant a la hauteur des thermiques.
-      do  ig=1,ngrid
-         zmax_sec(ig)=0.
-         zlevinter_dry(ig)=zlev(ig,1)
-      enddo
-      do  ig=1,ngrid
-! calcul de zlevinter
-          zlevinter_dry(ig)=zlev(ig,lmax_dry(ig)) + &
-     &    (linter_dry(ig)-lmax_dry(ig))*(zlev(ig,lmax_dry(ig)+1)-zlev(ig,lmax_dry(ig)))
-      zmax_sec(ig)=max(zmax_sec(ig),zlevinter_dry(ig)-zlev(ig,lmin(ig)))
-      enddo
-
-! ===========================================================================
-! ============= FIN DRY =====================================================
-! ===========================================================================
-
-
 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
 
@@ -857,5 +746,5 @@
 ! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
       llmax=1
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if (lalim(ig)>llmax) llmax=lalim(ig)
       enddo
@@ -865,7 +754,7 @@
 !   alim_star^2/(rho dz)
       do k=1,llmax-1
-         do ig=1,ngrid
+         do ig=1,ngridmx
             if (k<lalim(ig)) then
-               alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)**2  &
+         alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k)  &
       &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
          alim_star_tot_clos(ig)=alim_star_tot_clos(ig)+alim_star_clos(ig,k)
@@ -873,13 +762,16 @@
          enddo
       enddo
-
+ 
 ! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy
 ! True ratio is 3.5 but wetake into account the vmoy is the one alimentating
 ! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche)
 ! And r_aspect has been changed from 2 to 1.5 from observations
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if (alim_star2(ig)>1.e-10) then
-            f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/  &
-      &     (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
+!            f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/  &
+!      &     (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
+             f(ig)=wmax(ig)*alim_star_tot_clos(ig)/  &
+      &     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
+
          endif
       enddo
@@ -915,11 +807,11 @@
 !-------------------------------------------------------------------------
 
-      do l=1,nlay
+      do l=1,nlayermx
          entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
 
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
             if (l.lt.lmax(ig)) then
                fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
@@ -937,6 +829,6 @@
 ! autres corrections.
 
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
             if (detr(ig,l).gt.fm(ig,l)) then
                ncorecfm8=ncorecfm8+1
@@ -954,7 +846,7 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      do l=1,nlay
-
-         do ig=1,ngrid
+      do l=1,nlayermx
+
+         do ig=1,ngridmx
             if (l.lt.lmax(ig)) then
                fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
@@ -972,5 +864,5 @@
 !-------------------------------------------------------------------------
 
-         do ig=1,ngrid
+         do ig=1,ngridmx
             if (fm(ig,l+1).lt.0.) then
               print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
@@ -988,5 +880,5 @@
 !-------------------------------------------------------------------------
 !      if (iflag_thermals_optflux==0) then
-!         do ig=1,ngrid
+!         do ig=1,ngridmx
 !          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
 !     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
@@ -1007,5 +899,5 @@
 !-------------------------------------------------------------------------
 !      if (iflag_thermals_optflux==0) then
-!         do ig=1,ngrid
+!         do ig=1,ngridmx
 !            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
 !              f_old=fm(ig,l+1)
@@ -1021,5 +913,5 @@
 !-------------------------------------------------------------------------
 
-         do ig=1,ngrid
+         do ig=1,ngridmx
             if (detr(ig,l).gt.fm(ig,l)) then
                ncorecfm6=ncorecfm6+1
@@ -1042,5 +934,5 @@
 !-------------------------------------------------------------------------
 
-         do ig=1,ngrid
+         do ig=1,ngridmx
             if (fm(ig,l+1).lt.0.) then
                detr(ig,l)=detr(ig,l)+fm(ig,l+1)
@@ -1067,5 +959,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-        do ig=1,ngrid
+        do ig=1,ngridmx
            if (zw2(ig,l+1).gt.1.e-10) then
            zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
@@ -1088,6 +980,6 @@
 !-----------------------------------------------------------------------
 
-      do l=1,nlay-1
-         do ig=1,ngrid
+      do l=1,nlayermx-1
+         do ig=1,ngridmx
             eee0=entr(ig,l)
             ddd0=detr(ig,l)
@@ -1119,5 +1011,5 @@
 !              ddd=detr(ig)-entre
 !on s assure que tout s annule bien en zmax
-      do ig=1,ngrid
+      do ig=1,ngridmx
          fm(ig,lmax(ig)+1)=0.
          entr(ig,lmax(ig))=0.
@@ -1130,5 +1022,5 @@
 
 !IM 090508 beg
-      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then
+      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngridmx/4. ) then
          print*,'thermcell warning : large number of corrections'
          print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
@@ -1160,5 +1052,5 @@
       zdthladj(:,:)=0.
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if(lmax(ig) .gt. 1) then
          do k=1,lmax(ig)
@@ -1181,5 +1073,5 @@
       ztvd(:,:)=ztv(:,:)
       fm_down(:,:)=0.
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if (lmax(ig) .gt. 1) then
          do l=1,lmax(ig)
@@ -1219,5 +1111,5 @@
        zdthladj_down(:,:)=0.
 
-      do ig=1,ngrid
+      do ig=1,ngridmx
          if(lmax(ig) .gt. 1) then
          do k=1,lmax(ig)
@@ -1237,10 +1129,10 @@
 ! Calcul de la fraction de l'ascendance
 !------------------------------------------------------------------
-      do ig=1,ngrid
+      do ig=1,ngridmx
          fraca(ig,1)=0.
-         fraca(ig,nlay+1)=0.
-      enddo
-      do l=2,nlay
-         do ig=1,ngrid
+         fraca(ig,nlayermx+1)=0.
+      enddo
+      do l=2,nlayermx
+         do ig=1,ngridmx
             if (zw2(ig,l).gt.1.e-10) then
             fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
@@ -1273,6 +1165,6 @@
       nlarga0=0.
 
-      do k=1,nlay
-         do ig=1,ngrid
+      do k=1,nlayermx
+         do ig=1,ngridmx
             detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
          enddo
@@ -1280,5 +1172,5 @@
 
 !   calcul de la valeur dans les ascendances
-      do ig=1,ngrid
+      do ig=1,ngridmx
          zua(ig,1)=zu(ig,1)
          zva(ig,1)=zv(ig,1)
@@ -1287,7 +1179,7 @@
       enddo
 
-      gamma(1:ngrid,1)=0.
-      do k=2,nlay
-         do ig=1,ngrid
+      gamma(1:ngridmx,1)=0.
+      do k=2,nlayermx
+         do ig=1,ngridmx
             ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k)
             if(ltherm(ig,k).and.zmax(ig)>0.) then
@@ -1305,7 +1197,7 @@
       gamma(:,:)=0.
 
-      do k=2,nlay
-
-         do ig=1,ngrid
+      do k=2,nlayermx
+
+         do ig=1,ngridmx
 
             if (ltherm(ig,k)) then
@@ -1327,5 +1219,5 @@
 
       do iter=1,5
-         do ig=1,ngrid
+         do ig=1,ngridmx
 ! Pour memoire : calcul prenant en compte la fraction reelle
 !              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
@@ -1363,25 +1255,25 @@
       enddo
 
-      enddo ! k=2,nlay
+      enddo ! k=2,nlayermx
 
 ! Calcul du flux vertical de moment dans l'environnement.
 !---------------------------------------------------------
-      do k=2,nlay
-         do ig=1,ngrid
+      do k=2,nlayermx
+         do ig=1,ngridmx
             wud(ig,k)=fm(ig,k)*ue(ig,k)
             wvd(ig,k)=fm(ig,k)*ve(ig,k)
          enddo
       enddo
-      do ig=1,ngrid
+      do ig=1,ngridmx
          wud(ig,1)=0.
-         wud(ig,nlay+1)=0.
+         wud(ig,nlayermx+1)=0.
          wvd(ig,1)=0.
-         wvd(ig,nlay+1)=0.
+         wvd(ig,nlayermx+1)=0.
       enddo
 
 ! calcul des tendances.
 !-----------------------
-      do k=1,nlay
-         do ig=1,ngrid
+      do k=1,nlayermx
+         do ig=1,ngridmx
             pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k)  &
      &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
@@ -1412,8 +1304,8 @@
 
       modname='momentum'
-      call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
+      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
      &     masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax)
 
-      call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
+      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
      &     masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax)
 
@@ -1424,6 +1316,6 @@
 !------------------------------------------------------------------
 
-      do l=1,nlay
-         do ig=1,ngrid
+      do l=1,nlayermx
+         do ig=1,ngridmx
            pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
          enddo
@@ -1434,8 +1326,8 @@
 !------------------------------------------------------------------
 
-      if (nq .ne. 0.) then
+      if (nqmx .ne. 0.) then
       modname='tracer'
-      DO iq=1,nq
-          call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
+      DO iq=1,nqmx
+          call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
           &      masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax)
 
@@ -1448,5 +1340,5 @@
 
       modname='tke'
-      call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
+      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
       &      masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
 
@@ -1464,6 +1356,6 @@
 
     
-      do l=1,nlay-1
-       do ig=1,ngrid
+      do l=1,nlayermx-1
+       do ig=1,ngridmx
         teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))
         teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))
@@ -1471,14 +1363,14 @@
        enddo
       enddo
-      do ig=1,ngrid
-       teta_th_int(ig,nlay)=teta_th_int(ig,nlay-1)
-       teta_env_int(ig,nlay)=teta_env_int(ig,nlay-1)
-       teta_down_int(ig,nlay)=teta_down_int(ig,nlay-1)
-      enddo
-      do l=1,nlay
-       do ig=1,ngrid
+      do ig=1,ngridmx
+       teta_th_int(ig,nlayermx)=teta_th_int(ig,nlayermx-1)
+       teta_env_int(ig,nlayermx)=teta_env_int(ig,nlayermx-1)
+       teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
+      enddo
+      do l=1,nlayermx
+       do ig=1,ngridmx
         heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
-        buoyancyOut(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-        buoyancyEst(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
         heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
        enddo
