Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 627)
+++ /trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 628)
@@ -57,5 +57,5 @@
       REAL d_v_ajs(ngridmx,nlayermx) 
       REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
-      REAL detr_therm(ngridmx,nlayermx)
+      REAL detr_therm(ngridmx,nlayermx),detrmod(ngridmx,nlayermx)
       REAL zw2(ngridmx,nlayermx+1)
       REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
@@ -67,4 +67,5 @@
       LOGICAL qtransport_thermals,dtke_thermals
       INTEGER l,ig,iq,ii(1)
+      CHARACTER (LEN=20) modname
 
 !--------------------------------------------------------
@@ -87,4 +88,5 @@
       REAL zzw2(ngridmx,nlayermx+1)
       REAL zmax(ngridmx)
+      INTEGER ndt,zlmax
 
 !--------------------------------------------------------
@@ -225,20 +227,22 @@
 ! Initialization of intermediary variables
 
-         zfm_therm(:,:)=0.
-         zentr_therm(:,:)=0.
-         zdetr_therm(:,:)=0.
-         zheatFlux(:,:)=0.
-         zheatFlux_down(:,:)=0.
-         zbuoyancyOut(:,:)=0.
-         zbuoyancyEst(:,:)=0.
+!         zfm_therm(:,:)=0. !init is done inside
+!         zentr_therm(:,:)=0.
+!         zdetr_therm(:,:)=0.
+!         zheatFlux(:,:)=0.  
+!         zheatFlux_down(:,:)=0.
+!         zbuoyancyOut(:,:)=0.
+!         zbuoyancyEst(:,:)=0.
          zzw2(:,:)=0.
          zmax(:)=0.
          lmax(:)=0
-         d_t_the(:,:)=0.
-         d_u_the(:,:)=0.
-         d_v_the(:,:)=0.
+!         d_t_the(:,:)=0. !init is done inside
+
+!         d_u_the(:,:)=0. !transported outside
+!         d_v_the(:,:)=0.
          dq2_the(:,:)=0.
-         if (nqmx .ne. 0) then
-            d_q_the(:,:,:)=0.
+
+         if (nqmx .ne. 0 .and. ico2 .ne. 0) then
+            d_q_the(:,:,ico2)=0.
          endif
 
@@ -256,13 +260,13 @@
 
             d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact
-            d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
-            d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
-!            dq2_the(:,:)=dq2_the(:,:)*fact
+!            d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
+!            d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
+            dq2_the(:,:)=dq2_the(:,:)*fact
             if (ico2 .ne. 0) then
                d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact
             endif
 
-             zmaxth(:)=zmaxth(:)+zmax(:)*fact
-             lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
+            zmaxth(:)=zmaxth(:)+zmax(:)*fact
+            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
             fm_therm(:,:)=fm_therm(:,:)  &
      &      +zfm_therm(:,:)*fact
@@ -288,6 +292,6 @@
 
            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
-           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
-           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
+!           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
+!           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
             if (ico2 .ne. 0) then
                d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
@@ -297,6 +301,6 @@
 
             zt(:,:) = zt(:,:) + d_t_the(:,:)
-            zu(:,:) = zu(:,:) + d_u_the(:,:)
-            zv(:,:) = zv(:,:) + d_v_the(:,:)
+!            zu(:,:) = zu(:,:) + d_u_the(:,:)
+!            zv(:,:) = zv(:,:) + d_v_the(:,:)
             if (ico2 .ne. 0) then
              pq_therm(:,:,ico2) = &
@@ -310,4 +314,10 @@
 
       lmax(:)=nint(lmax_real(:))
+      zlmax=MAXVAL(lmax(:))+2
+      if (zlmax .ge. nlayermx) then
+        print*,'thermals have reached last layer of the model'
+        print*,'this is not good !'
+      endif
+
 
 ! Now that we have computed total entrainment and detrainment, we can
@@ -321,11 +331,23 @@
       enddo
 
-!      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
-!     &      ,fm_therm,entr_therm,  &
-!     &     masse,zu,d_u_ajs)
-!
-!      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
-!     &       ,fm_therm,entr_therm,  &
-!     &     masse,zv,d_v_ajs)
+      detrmod(:,:)=0.
+      do l=1,zlmax
+         do ig=1,ngridmx
+            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
+     &      +entr_therm(ig,l)
+            if (detrmod(ig,l).lt.0.) then
+               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
+               detrmod(ig,l)=0.
+            endif
+         enddo
+      enddo
+      ndt=10
+      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
+     &      ,fm_therm,entr_therm,detrmod,  &
+     &     masse,zu,d_u_ajs,ndt,zlmax)
+
+      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
+     &       ,fm_therm,entr_therm,detrmod,  &
+     &     masse,zv,d_v_ajs,ndt,zlmax)
 
       if (nqmx .ne. 0.) then
@@ -333,6 +355,6 @@
       if (iq .ne. ico2) then
       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
-     &     ,fm_therm,entr_therm,  &
-     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq))
+     &     ,fm_therm,entr_therm,detrmod,  &
+     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,zlmax)
       endif
       ENDDO
@@ -340,7 +362,19 @@
 
       if (dtke_thermals) then
+      detrmod(:,:)=0.
+      ndt=1
+      do l=1,zlmax
+         do ig=1,ngridmx
+            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
+     &      +entr_therm(ig,l)
+            if (detrmod(ig,l).lt.0.) then
+               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
+               detrmod(ig,l)=0.
+            endif
+         enddo
+      enddo
       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
-     &     ,fm_therm,entr_therm,  &
-     &    masse,q2_therm,dq2_therm)
+     &     ,fm_therm,entr_therm,detrmod,  &
+     &    masse,q2_therm,dq2_therm,ndt,zlmax)
       endif
 
@@ -363,6 +397,8 @@
 ! **********************************************************************
 
-           pdu_th(:,:)=d_u_ajs(:,:)/ptimestep
-           pdv_th(:,:)=d_v_ajs(:,:)/ptimestep
+      do l=1,zlmax
+        pdu_th(:,l)=d_u_ajs(:,l)
+        pdv_th(:,l)=d_v_ajs(:,l)
+      enddo
 
            if(qtransport_thermals) then
@@ -370,7 +406,11 @@
                do iq=1,nqmx
                 if (iq .ne. ico2) then
-                  pdq_th(:,:,iq)=d_q_ajs(:,:,iq)
+                  do l=1,zlmax
+                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)
+                  enddo
                 else
-                  pdq_th(:,:,iq)=d_q_ajs(:,:,iq)/ptimestep
+                  do l=1,zlmax
+                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep
+                  enddo
                 endif
                enddo
@@ -387,5 +427,7 @@
            ENDIF
 
-           pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
+           do l=1,zlmax
+              pdt_th(:,l)=d_t_ajs(:,l)/ptimestep
+           enddo
 
 
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90	(revision 627)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90	(revision 628)
@@ -1,7 +1,6 @@
-      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
-     &    masse0,q_therm,dq_therm)
+      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm,entr,detr,  &
+     &    masse0,q_therm,dq_therm,ndt,zlmax)
       implicit none
 
-! #include "iniprint.h"
 !=======================================================================
 !
@@ -11,4 +10,7 @@
 !   Version modifiee pour prendre les downdrafts a la place de la 
 !   subsidence compensatoire
+!
+!   Version with sub-timestep for Martian thin layers
+!
 !=======================================================================
 
@@ -20,8 +22,11 @@
       INTEGER, INTENT(IN) :: ngrid,nlayer
       REAL, INTENT(IN) :: ptimestep
-      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
-      REAL, INTENT(IN) ::entr0(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: fm(ngridmx,nlayermx+1)
+      REAL, INTENT(IN) :: entr(ngridmx,nlayermx)
+      REAL, INTENT(IN) :: detr(ngridmx,nlayermx)
       REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
       REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
+      INTEGER, INTENT(IN) :: ndt
+      INTEGER, INTENT(IN) :: zlmax
 
 ! ============================ OUTPUTS ===========================
@@ -33,8 +38,7 @@
       REAL q(ngridmx,nlayermx)
       REAL qa(ngridmx,nlayermx)
-      INTEGER ig,k
-      REAL detr(ngridmx,nlayermx),entr(ngridmx,nlayermx)
-      REAL wqd(ngridmx,nlayermx+1)
-      REAL zzm
+      INTEGER ig,k,i
+      REAL invflux0(ngridmx,nlayermx)
+      REAL ztimestep
 
 ! =========== Init ==============================================
@@ -42,18 +46,12 @@
       qa(:,:)=q_therm(:,:)
       q(:,:)=q_therm(:,:)
-      entr(:,:)=entr0(:,:)
-      detr(:,:)=0.
 
-      do k=1,nlayermx
-         do ig=1,ngridmx
-            detr(ig,k)=fm0(ig,k)-fm0(ig,k+1)+entr(ig,k)
-            if (detr(ig,k).lt.0.) then
-               entr(ig,k)=entr(ig,k)-detr(ig,k)
-               detr(ig,k)=0.
-            endif
-         enddo
-      enddo
+! ====== Computing q ============================================
 
-! =========== Updraft ============================================
+      dq_therm(:,:)=0.
+      ztimestep=ptimestep/real(ndt)
+      invflux0(:,:)=ztimestep/masse0(:,:)      
+
+      do i=1,ndt
 
       do ig=1,ngrid
@@ -61,10 +59,10 @@
       enddo
 
-      do k=2,nlayermx
+      do k=2,zlmax
          do ig=1,ngridmx
-            if ((fm0(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
+            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
      &         1.e-5*masse0(ig,k)) then
-         qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
-     &     /(fm0(ig,k+1)+detr(ig,k))
+         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
+     &     /(fm(ig,k+1)+detr(ig,k))
             else
                qa(ig,k)=q(ig,k)
@@ -73,39 +71,23 @@
       enddo
 
-! =========== Subsidence ==========================================
-
-      do k=2,nlayermx
-         do ig=1,ngridmx
-
-! Schema avec advection sur plus qu'une maille.
-            zzm=masse0(ig,k)/ptimestep
-            if (fm0(ig,k)>zzm) then
-               wqd(ig,k)=(zzm*q(ig,k)+(fm0(ig,k)-zzm)*q(ig,k+1))
-            else
-               wqd(ig,k)=fm0(ig,k)*q(ig,k)
-            endif
-         enddo
-      enddo
-      do ig=1,ngrid
-         wqd(ig,1)=0.
-         wqd(ig,nlayermx+1)=0.
+      do k=1,zlmax
+           q(:,k)=q(:,k)+         &
+     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
+     &    -fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1))  &
+     &               *invflux0(:,k)
       enddo
 
+      enddo !of do i=1,ndt
 
-! ====== dq ======================================================
+! ====== Derivative ==============================================
 
-      dq_therm(:,:)=0.
 
-         do k=1,nlayermx-1
-           q(:,k)=q(:,k)+         &
-     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
-     &       -wqd(:,k)+wqd(:,k+1))  &
-     &               *ptimestep/masse0(:,k)
-         enddo
-
-         do k=1,nlayermx-1
+         do k=1,zlmax
           dq_therm(:,k)=(q(:,k)-q_therm(:,k))/ptimestep
          enddo
 
+! ==============
+
       return
       end
+
Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 627)
+++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 628)
@@ -99,10 +99,12 @@
       REAL f(ngridmx)
 
+      REAL detrmod(ngridmx,nlayermx)
+
       REAL teta_th_int(ngridmx,nlayermx)
       REAL teta_env_int(ngridmx,nlayermx)
       REAL teta_down_int(ngridmx,nlayermx)
 
-      CHARACTER (LEN=20) :: modname
       CHARACTER (LEN=80) :: abort_message
+      INTEGER ndt
 
 ! ============= PLUME VARIABLES ============
@@ -123,4 +125,5 @@
       REAL denom(ngridmx)
       REAL zlevinter(ngridmx)
+      INTEGER zlmax
 
 ! =========================================
@@ -186,7 +189,8 @@
       detr(:,:)=0.
       fm(:,:)=0.
-      zu(:,:)=pu(:,:)
-      zv(:,:)=pv(:,:)
+!      zu(:,:)=pu(:,:)
+!      zv(:,:)=pv(:,:)
       zhc(:,:)=pt(:,:)/zpopsk(:,:)
+      ndt=1
 
 ! **********************************************************************
@@ -888,4 +892,10 @@
 ! ===========================================================================
 
+      zlmax=MAXVAL(lmax(:))+2
+      if (zlmax .ge. nlayermx) then
+        print*,'thermals have reached last layer of the model'
+        print*,'this is not good !'
+      endif
+
 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
 
@@ -968,10 +978,10 @@
 !-------------------------------------------------------------------------
 
-      do l=1,nlayermx
+      do l=1,zlmax
          entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
 
-      do l=1,nlayermx
+      do l=1,zlmax
          do ig=1,ngridmx
             if (l.lt.lmax(ig)) then
@@ -1007,5 +1017,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      do l=1,nlayermx
+      do l=1,zlmax
 
          do ig=1,ngridmx
@@ -1156,5 +1166,5 @@
 !-----------------------------------------------------------------------
 
-      do l=1,nlayermx-1
+      do l=1,zlmax
          do ig=1,ngridmx
             eee0=entr(ig,l)
@@ -1347,9 +1357,6 @@
 ! Calcul de la fraction de l'ascendance
 !------------------------------------------------------------------
-      do ig=1,ngridmx
-         fraca(ig,1)=0.
-         fraca(ig,nlayermx+1)=0.
-      enddo
-      do l=2,nlayermx
+      fraca(:,:)=0.
+      do l=2,zlmax
          do ig=1,ngridmx
             if (zw2(ig,l).gt.1.e-10) then
@@ -1520,12 +1527,24 @@
 
       else
-
-      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
-     &      ,fm,entr,  &
-     &     masse,zu,pduadj)
-
-      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
-     &       ,fm,entr,  &
-     &     masse,zv,pdvadj)
+!      detrmod(:,:)=0.
+!      do k=1,zlmax
+!         do ig=1,ngridmx
+!            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
+!     &      +entr(ig,k)
+!            if (detrmod(ig,k).lt.0.) then
+!               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
+!               detrmod(ig,k)=0.
+!            endif
+!         enddo
+!      enddo
+!
+!
+!      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
+!     &      ,fm,entr,detrmod,  &
+!     &     masse,zu,pduadj,ndt,zlmax)
+!
+!      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
+!     &       ,fm,entr,detrmod,  &
+!     &     masse,zv,pdvadj,ndt,zlmax)
 
       endif
@@ -1538,22 +1557,32 @@
 ! The rest is transported outside the sub-timestep loop
 
+      ratiom(:,:)=1.
+
       if (ico2.ne.0) then
-!      if (nqmx .ne. 0.) then
+      detrmod(:,:)=0.
+      do k=1,zlmax
+         do ig=1,ngridmx
+            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
+     &      +entr(ig,k)
+            if (detrmod(ig,k).lt.0.) then
+               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
+               detrmod(ig,k)=0.
+            endif
+         enddo
+      enddo
+
       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
-     &     ,fm,entr,  &
-     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2))
-!      endif
+     &     ,fm,entr,detrmod,  &
+     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),ndt,zlmax)
 
 ! Compute the ratio between theta and theta_m
       
-       do l=1,nlayermx
+       do l=1,zlmax
           do ig=1,ngridmx
              ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
           enddo
        enddo
-       else
-          ratiom(:,:)=1.
+
        endif
-
 
 !------------------------------------------------------------------
@@ -1561,5 +1590,6 @@
 !------------------------------------------------------------------
 
-      do l=1,nlayermx
+      pdtadj(:,:)=0.
+      do l=1,zlmax
          do ig=1,ngridmx
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
@@ -1601,5 +1631,9 @@
        teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
       enddo
-      do l=1,nlayermx
+        heatFlux(:,:)=0.
+        buoyancyOut(:,:)=0.
+        buoyancyEst(:,:)=0.
+        heatFlux_down(:,:)=0.
+      do l=1,zlmax
        do ig=1,ngridmx
         heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
