Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1020)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1028)
@@ -12,5 +12,5 @@
      &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
      &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
-     &   ,tifeedback,microphys,caps,photochem,calltherm,outptherm       &
+     &   ,tifeedback,microphys,caps,photochem,calltherm                 &
      &   ,callrichsl,callslope,tituscap
      
@@ -26,5 +26,5 @@
      &   ,callnirco2,callnlte,callthermos,callconduct,                  &
      &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
-     &   ,calltherm,outptherm,callrichsl,callslope,tituscap
+     &   ,calltherm,callrichsl,callslope,tituscap
 
 
Index: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 1020)
+++ trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 1028)
@@ -1,5 +1,18 @@
+! ---------------------------------------------------------------------
+! Arnaud Colaitis 2011-01-05
 !
-! AC 2011-01-05
+! Interface routine between physiq.F (driver for physics) and thermcell_main_mars (thermals model)
+! Handles sub-timestep for the thermals model, outputs and diagnostics
+! 
+! The thermals model implemented in the Mars LMD-GCM is called after the vertical turbulent mixing scheme (Mellor and Yamada)  
+! and surface layer scheme (Richardson-based surface layer with subgrid gustiness parameterization)
+! Other routines called before the thermals model are Radiative transfer scheme (Madeleine et.al)  
+! and gravity wave and subgrid scale topography drag.
 !
+! Reference paper for the Martian thermals model is Colaitis et al. 2013 (JGR-planets)
+! Mentions in the routines to a "paper" points to the above reference.
+!
+! This model is based on the LMDZ thermals model from C. Rio and F. Hourdin
+! ---------------------------------------------------------------------
       SUBROUTINE calltherm_interface (firstcall, & 
      & zzlev,zzlay, &
@@ -9,12 +22,21 @@
      & pdhdif,hfmax,wstar,sensibFlux)
 
-       USE ioipsl_getincom
-
+
+      USE ioipsl_getincom, only : getin ! to read options from external file "run.def"
       implicit none
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "tracer.h"
+#include "callkeys.h" !contains logical values determining which subroutines and which options are used. 
+                      ! includes logical "tracer", which is .true. if tracers are present and to be transported
+                      ! includes integer "iradia", frequency of calls to radiative scheme wrt calls to physics (and is typically==1)
+#include "dimensions.h" !contains global GCM grid dimension informations
+#include "dimphys.h" !similar to dimensions.h, and based on it; includes
+                     ! "ngridmx": number of horizontal grid points
+                     ! "nlayermx": number of atmospheric layers
+#include "comcstfi.h" !contains physical constant values such as
+                      ! "g" : gravitational acceleration (m.s-2)
+                      ! "r" : recuced gas constant (J.K-1.mol-1)
+                      ! "cpp" : specific heat of the atmosphere (J.kg-1.K-1)
+#include "tracer.h" !contains tracer information such as
+                    ! "nqmx": number of tracers
+                    ! "noms()": tracer name
 
 !--------------------------------------------------------
@@ -22,21 +44,24 @@
 !--------------------------------------------------------
 
-!      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
-      REAL, INTENT(IN) :: ptimestep
-      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
-      REAL, INTENT(IN) :: 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(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
-      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx)
-      REAL, INTENT(IN) :: pdt(ngridmx,nlayermx)
-      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
-      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
-      REAL, INTENT(IN) :: pdhdif(ngridmx,nlayermx)
-      REAL, INTENT(IN) :: sensibFlux(ngridmx)
+      REAL, INTENT(IN) :: ptimestep !timestep (s)
+      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) !intermediate pressure levels (Pa)
+      REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) !Pressure at the middle of the layers (Pa)
+      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) !Geopotential at the middle of the layers (m2s-2)
+      REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx) !u,v components of the wind (ms-1)
+      REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)!temperature (K) and tracer concentration (kg/kg)
+      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
+      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      LOGICAL, INTENT(IN) :: firstcall ! =.true. if this is the first time the routine is called. 
+                                       ! Used to detect and store which tracer is co2
+      REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx) ! wind velocity change from routines called 
+                                                                      ! before thermals du/dt (m/s/s)
+      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx) ! tracer concentration change from routines called 
+                                                     ! before thermals dq/dt (kg/kg/s)
+      REAL, INTENT(IN) :: pdt(ngridmx,nlayermx) ! temperature change from routines called before thermals dT/dt (K/s)
+      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1) ! turbulent kinetic energy
+      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx) ! ratio of pressure at middle of layer to surface pressure, 
+                                                   ! to the power r/cp, i.e. zpopsk=(pplay(ig,l)/pplev(ig,1))**rcp
+      REAL, INTENT(IN) :: pdhdif(ngridmx,nlayermx) ! potential temperature change from turbulent diffusion scheme dT/dt (K/s)
+      REAL, INTENT(IN) :: sensibFlux(ngridmx) ! sensible heat flux computed from surface layer scheme
 
 !--------------------------------------------------------
@@ -44,12 +69,21 @@
 !--------------------------------------------------------
 
-      REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx)
-      REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx)
-      REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx)
-      REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx)
-      INTEGER, INTENT(OUT) :: lmax(ngridmx)
-      REAL, INTENT(OUT) :: zmaxth(ngridmx)
-      REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1)
-      REAL, INTENT(OUT) :: wstar(ngridmx)
+      REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx) ! wind velocity change from thermals du/dt (m/s/s)
+      REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx) ! wind velocity change from thermals dv/dt (m/s/s)
+      REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx) ! temperature change from thermals dT/dt (K/s)
+      REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx) ! tracer change from thermals dq/dt (kg/kg/s)
+      INTEGER, INTENT(OUT) :: lmax(ngridmx) ! layer number reacher by thermals in grid point
+      REAL, INTENT(OUT) :: zmaxth(ngridmx) ! equivalent to lmax, but in (m), interpolated
+      REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1) ! turbulent kinetic energy change from thermals dtke/dt
+      REAL, INTENT(OUT) :: wstar(ngridmx) ! free convection velocity (m/s)
+
+!--------------------------------------------------------
+! Control parameters of the thermal model
+!--------------------------------------------------------
+! such variables are: 
+!  - outptherm to control internal diagnostics
+!  - qtransport_thermals to control tracer transport in thermals
+!  - dtke_thermals to control turbulent kinetic energy transport in thermals
+     LOGICAL,SAVE :: outptherm,qtransport_thermals,dtke_thermals
 
 !--------------------------------------------------------
@@ -65,10 +99,9 @@
       REAL zw2(ngridmx,nlayermx+1)
       REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
-      REAL ztla(ngridmx,nlayermx)
       REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
       REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
       REAL lmax_real(ngridmx)
       REAL masse(ngridmx,nlayermx)
-      LOGICAL qtransport_thermals,dtke_thermals
+
       INTEGER l,ig,iq,ii(1),k
       CHARACTER (LEN=20) modname
@@ -79,6 +112,4 @@
 
       REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
-      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
-      REAL dq2_the(ngridmx,nlayermx)
       INTEGER isplit
       INTEGER,SAVE :: nsplit_thermals
@@ -113,9 +144,37 @@
 !--------------------------------------------------------
 
-      INTEGER ico2
-      SAVE ico2
-
-! **********************************************************************
-! Initialization
+      INTEGER,SAVE :: ico2
+
+
+      if(firstcall) then
+        ico2=0
+! **********************************************************************
+! Polar night mixing : theta_m
+! Get index of co2 tracer in tracer array
+! **********************************************************************
+        if (tracer) then !NB: tracer==.true. if tracers are present and to be transported
+!     Prepare Special treatment if one of the tracers is CO2 gas
+           do iq=1,nqmx
+             if (noms(iq).eq."co2") then
+                ico2=iq
+             end if
+           enddo
+        endif
+      
+        ! control parameters of the thermal model:
+        qtransport_thermals=.true. !! default setting: transport tracers in thermals
+        
+        dtke_thermals=.false. !! default setting
+        !  switch to transport TKE in thermals. still experimental, for testing purposes only. not used in current thermal plume models both on Earth and Mars.
+
+        outptherm=.false. !! default setting
+        ! switch to control if internal quantities ofthe thermals must be output (typically set to .false., but very useful for fine diagnostics. 
+        ! get value (if set by user) of this parameter from run.def file
+        call getin("outptherm",outptherm)
+
+      endif !of if firstcall
+
+! **********************************************************************
+! Initializations
 ! **********************************************************************
 
@@ -128,5 +187,4 @@
       q2_therm(:,:)=0.
       dq2_therm(:,:)=0.
-      ztla(:,:)=0.
       pbl_dtke(:,:)=0.
       fm_therm(:,:)=0.
@@ -150,23 +208,21 @@
 
 ! **********************************************************************
-! Preparing inputs for the thermals
-! **********************************************************************
-
-       zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
-       zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
-       zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
-
-       pq_therm(:,:,:)=0.
-       qtransport_thermals=.true. !! default setting
-       !call getin("qtransport_thermals",qtransport_thermals)
+! Preparing inputs for the thermals: increment tendancies 
+! from other subroutines called before the thermals model
+! **********************************************************************
+
+       zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep ! u-component of wind velocity
+       zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep ! v-component of wind velocity
+       zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep ! temperature
+
+       pq_therm(:,:,:)=0. ! tracer concentration
 
        if(qtransport_thermals) then
-          if(tracer) then
-                pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
+          if(tracer) then ! tracer is a logical that is true if tracer must be transported in the GCM physics
+                pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep ! tracer concentration
           endif
        endif
 
-       dtke_thermals=.false. !! default setting
-       call getin("dtke_thermals",dtke_thermals)
+
        IF(dtke_thermals) THEN
           DO l=1,nlayermx
@@ -174,20 +230,4 @@
           ENDDO
        ENDIF
-
-! **********************************************************************
-! Polar night mixing : theta_m
-! **********************************************************************
-
-      if(firstcall) then
-        ico2=0
-        if (tracer) then
-!     Prepare Special treatment if one of the tracers is CO2 gas
-           do iq=1,nqmx
-             if (noms(iq).eq."co2") then
-                ico2=iq
-             end if
-           enddo
-        endif
-      endif !of if firstcall
 
 
@@ -206,20 +246,22 @@
                                 ! chosen timestep for the radiative transfer.
                                 ! It is recommended to run with 96 timestep per day and 
-                                ! iradia = 1., configuration in which thermals can run
+                                ! call radiative transfer at each timestep, 
+                                ! configuration in which thermals can run
                                 ! very well with a sub-timestep of 10.
          IF (firstcall) THEN
             r_aspect_thermals=1.  ! same value is OK for GCM and mesoscale
+! Sub Timestep for the mesoscale model:
 #ifdef MESOSCALE
             !! valid for timesteps < 200s
             nsplit_thermals=4
 #else
+! Sub Timestep for the GCM:
+            ! If there is at least 96 timesteps per day in the gcm, subtimestep factor 10
             IF ((ptimestep .le. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
                nsplit_thermals=10
-            ELSE
+            ELSE !if there is less, subtimestep factor 35
                nsplit_thermals=35
             ENDIF
 #endif
-            call getin("nsplit_thermals",nsplit_thermals)
-            call getin("r_aspect_thermals",r_aspect_thermals)
          ENDIF
 
@@ -228,95 +270,73 @@
 ! **********************************************************************
 
-         zdt=ptimestep/REAL(nsplit_thermals)
-
-         DO isplit=1,nsplit_thermals
+         zdt=ptimestep/REAL(nsplit_thermals) !subtimestep
+
+         DO isplit=1,nsplit_thermals !temporal loop on the subtimestep
 
 ! Initialization of intermediary variables
 
-!         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. !init is done inside
-
-!         d_u_the(:,:)=0. !transported outside
-!         d_v_the(:,:)=0.
-         dq2_the(:,:)=0.
-
-         if (nqmx .ne. 0 .and. ico2 .ne. 0) then
+
+         if (nqmx .ne. 0 .and. ico2 .ne. 0) then !initialize co2 tracer tendancy
             d_q_the(:,:,ico2)=0.
          endif
 
+! CALL to main thermal routine
              CALL thermcell_main_mars(zdt  &
      &      ,pplay,pplev,pphi,zzlev,zzlay  &
      &      ,zu,zv,zt,pq_therm,q2_therm  &
-     &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
+     &      ,d_t_the,d_q_the  &
      &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
      &      ,r_aspect_thermals &
      &      ,zzw2,fraca,zpopsk &
-     &      ,ztla,zheatFlux,zheatFlux_down &
+     &      ,zheatFlux,zheatFlux_down &
      &      ,zbuoyancyOut,zbuoyancyEst)
 
       fact=1./REAL(nsplit_thermals)
 
-            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
+! Update thermals tendancies
+            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact !temperature
             if (ico2 .ne. 0) then
-               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact
+               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact !co2 mass mixing ratio
             endif
-
-            zmaxth(:)=zmaxth(:)+zmax(:)*fact
-            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
-            fm_therm(:,:)=fm_therm(:,:)  &
+            zmaxth(:)=zmaxth(:)+zmax(:)*fact !thermals height
+            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact !thermals height index
+            fm_therm(:,:)=fm_therm(:,:)  & !upward mass flux
      &      +zfm_therm(:,:)*fact
-            entr_therm(:,:)=entr_therm(:,:)  &
+            entr_therm(:,:)=entr_therm(:,:)  & !entrainment mass flux
      &       +zentr_therm(:,:)*fact
-            detr_therm(:,:)=detr_therm(:,:)  &
+            detr_therm(:,:)=detr_therm(:,:)  & !detrainment mass flux
      &       +zdetr_therm(:,:)*fact
-            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
-
-            heatFlux(:,:)=heatFlux(:,:) &
+            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact !updraft fractional coverage
+            heatFlux(:,:)=heatFlux(:,:) & !upward heat flux
      &       +zheatFlux(:,:)*fact
-            heatFlux_down(:,:)=heatFlux_down(:,:) &
+            heatFlux_down(:,:)=heatFlux_down(:,:) & !downward heat flux
      &       +zheatFlux_down(:,:)*fact
-            buoyancyOut(:,:)=buoyancyOut(:,:) &
+            buoyancyOut(:,:)=buoyancyOut(:,:) & !plume final buoyancy
      &       +zbuoyancyOut(:,:)*fact
-            buoyancyEst(:,:)=buoyancyEst(:,:) &
+            buoyancyEst(:,:)=buoyancyEst(:,:) & !plume estimated buoyancy used for vertical velocity computation
      &       +zbuoyancyEst(:,:)*fact
-  
-
-            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
-
-!  accumulation de la tendance
-
-           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(:,:)
+            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact !vertical velocity
+
+! Save tendancies
+           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) !temperature tendancy (delta T)
             if (ico2 .ne. 0) then
-               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
+               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2) !tracer tendancy (delta q)
             endif
-!            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
-!  incrementation des variables meteo
-
-            zt(:,:) = zt(:,:) + d_t_the(:,:)
-!            zu(:,:) = zu(:,:) + d_u_the(:,:)
-!            zv(:,:) = zv(:,:) + d_v_the(:,:)
+
+!  Increment temperature and co2 concentration for next pass in subtimestep loop
+            zt(:,:) = zt(:,:) + d_t_the(:,:) !temperature
             if (ico2 .ne. 0) then
              pq_therm(:,:,ico2) = &
-     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2)
+     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2) !co2 tracer
             endif
-!            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
 
 
          ENDDO ! isplit
 !****************************************************************
+
+!     Do we have thermals that are too high ? 
 
       lmax(:)=nint(lmax_real(:))
@@ -329,6 +349,8 @@
 
 ! Now that we have computed total entrainment and detrainment, we can
-! advect u, v, and q in thermals. (theta already advected). We can do
-! that separatly because u,v,and q are not used in thermcell_main for
+! advect u, v, and q in thermals. (potential temperature and co2 MMR 
+! have already been advected in thermcell_main because they are coupled
+! to the determination of the thermals caracteristics). This is done 
+! separatly because u,v, and q are not used in thermcell_main for
 ! any thermals-related computation : they are purely passive.
 
@@ -338,4 +360,6 @@
       enddo
 
+! recompute detrainment mass flux from entrainment and upward mass flux
+! this ensure mass flux conservation
       detrmod(:,:)=0.
       do l=1,zlmax
@@ -349,4 +373,7 @@
          enddo
       enddo
+
+      ! u component of wind velocity advection in thermals
+      ! result is a derivative (d_u_ajs in m/s/s)
       ndt=10
       call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
@@ -354,7 +381,12 @@
      &     masse,zu,d_u_ajs,ndt,zlmax)
 
+      ! v component of wind velocity advection in thermals
+      ! result is a derivative (d_v_ajs in m/s/s)
       call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
      &       ,fm_therm,entr_therm,detrmod,  &
      &     masse,zv,d_v_ajs,ndt,zlmax)
+
+      ! non co2 tracers advection in thermals
+      ! result is a derivative (d_q_ajs in kg/kg/s)
 
       if (nqmx .ne. 0.) then
@@ -368,17 +400,7 @@
       endif
 
+      ! tke advection in thermals
+      ! result is a tendancy (d_u_ajs in J)
       if (dtke_thermals) then
-      detrmod(:,:)=0.
-      ndt=10
-      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,detrmod,  &
@@ -386,4 +408,5 @@
       endif
 
+      ! compute wmax for diagnostics
       DO ig=1,ngridmx
          wmax(ig)=MAXVAL(zw2(ig,:))
@@ -408,4 +431,5 @@
       enddo
 
+           ! if tracers are transported in thermals, update output variables, else these are 0.
            if(qtransport_thermals) then
               if(tracer) then
@@ -413,9 +437,9 @@
                 if (iq .ne. ico2) then
                   do l=1,zlmax
-                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)
+                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq) !non-co2 tracers d_q_ajs are dq/dt (kg/kg/s)
                   enddo
                 else
                   do l=1,zlmax
-                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep
+                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep !co2 tracer d_q_ajs is dq (kg/kg)
                   enddo
                 endif
@@ -424,4 +448,5 @@
            endif
 
+           ! if tke is transported in thermals, update output variable, else this is 0.
            IF(dtke_thermals) THEN
               DO l=2,nlayermx
@@ -433,4 +458,5 @@
            ENDIF
 
+           ! update output variable for temperature. d_t_ajs is delta T in (K), pdt_th is dT/dt in (K/s)
            do l=1,zlmax
               pdt_th(:,l)=d_t_ajs(:,l)/ptimestep
@@ -439,7 +465,11 @@
 
 ! **********************************************************************
-! Compute the free convection velocity scale for vdifc
-! **********************************************************************
-
+! SURFACE LAYER INTERFACE
+! Compute the free convection velocity w* scale for surface layer gustiness
+! speed parameterization. The computed value of w* will be used at the next
+! timestep to modify surface-atmosphere exchange fluxes, because the surface
+! layer scheme and diffusion are called BEFORE the thermals. (outside of these
+! routines)
+! **********************************************************************
 
 ! Potential temperature gradient
@@ -453,5 +483,5 @@
       ENDDO
 
-! Computation of the pbl mixed layer temperature
+! Computation of the PBL mixed layer temperature
 
       DO ig=1, ngridmx
@@ -460,5 +490,8 @@
       ENDDO
 
-! we must add the heat flux from the diffusion scheme to hfmax
+! In order to have an accurate w*, we must add the heat flux from the
+! diffusion scheme to the computation of the maximum heat flux hfmax
+! Here pdhdif is the potential temperature change from the diffusion
+! scheme (Mellor and Yamada, see paper section 6, paragraph 57) 
 
 ! compute rho as it is after the diffusion
@@ -482,11 +515,14 @@
       ENDDO
 
-! compute w'teta' from diffusion
+! compute w'theta' (vertical turbulent flux of temperature) from 
+! the diffusion scheme
 
       wtdif(:,:)=rpdhd(:,:)/rho(:,:)
 
+! Now we compute the contribution of the thermals to w'theta':
 ! compute rho as it is after the thermals
 
       rho(:,:)=pplay(:,:)/(r*(zt(:,:)))
+
 ! integrate -rho*pdhdif
 
@@ -508,9 +544,12 @@
       wtth(:,:)=rpdhd(:,:)/rho(:,:)
 
-! We get the max heat flux from thermals and add the contribution from the diffusion
+! Add vertical turbulent heat fluxes from the thermals and the diffusion scheme
+! and compute the maximum
 
       DO ig=1,ngridmx
         hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:))
       ENDDO
+
+! Finally we can compute the free convection velocity scale
 ! We follow Spiga et. al 2010 (QJRMS)
 ! ------------
@@ -524,12 +563,14 @@
       ENDDO
 
-
-
 ! **********************************************************************
 ! Diagnostics
 ! **********************************************************************
+! THESE DIAGNOSTICS ARE WRITTEN IN THE LMD-GCM FORMAT
+! THESE ARE LEFT AS A FURTHER INDICATION OF THE QUANTITIES DEFINED IN THIS ROUTINE
         
-        if(outptherm) then
-        if (ngridmx .eq. 1) then
+! outptherm is a logical in callkeys that controls if internal quantities of the thermals must be output.
+
+        if (outptherm) then
+        if (ngridmx .eq. 1) then !these are 1D outputs
         call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
      &                       'kg/m-2',1,entr_therm)
@@ -565,5 +606,5 @@
      &         'heat flux PBL','K.m/s',1,wtdif(:,:)+wtth(:,:))
 
-      else
+      else !these are 3D outputs
 
         call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
@@ -583,5 +624,5 @@
 
       endif
-      endif
+      endif ! of if (outptherm)
 
        END
Index: trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 1020)
+++ trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 1028)
@@ -267,9 +267,4 @@
          write(*,*) " calltherm = ",calltherm
 
-         write(*,*) "output thermal diagnostics ?"
-         outptherm=.false. ! default value
-         call getin("outptherm",outptherm)
-         write(*,*) " outptherm = ",outptherm
-
          write(*,*) "call convective adjustment ?"
          calladj=.true. ! default value
Index: trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90	(revision 1020)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90	(revision 1028)
@@ -1,2 +1,5 @@
+! ---------------------------------------------------------------------
+! Arnaud Colaitis 2011-01-05
+! --------------------------------------------------------------------
       subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm,entr,detr,  &
      &    masse0,q_therm,dq_therm,ndt,zlmax)
@@ -5,9 +8,8 @@
 !=======================================================================
 !
-!   Calcul du transport verticale dans la couche limite en presence
-!   de "thermiques" explicitement representes
-!   calcul du dq/dt une fois qu'on connait les ascendances
-!   Version modifiee pour prendre les downdrafts a la place de la 
-!   subsidence compensatoire
+!   Compute the thermals contribution of explicit thermals 
+!   to vertical transport in the PBL.
+!   dq is computed once upward, entrainment and detrainment mass fluxes
+!   are known.
 !
 !   Version with sub-timestep for Martian thin layers
@@ -15,18 +17,20 @@
 !=======================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+#include "dimensions.h" !contains global GCM grid dimension informations
+#include "dimphys.h" !similar to dimensions.h, and based on it; includes
+                     ! "ngridmx": number of horizontal grid points
+                     ! "nlayermx": number of atmospheric layers
 
 ! ============================ INPUTS ============================
 
-      INTEGER, INTENT(IN) :: ngrid,nlayer
-      REAL, INTENT(IN) :: ptimestep
-      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
+      INTEGER, INTENT(IN) :: ngrid,nlayer ! number of grid points and number of levels
+      REAL, INTENT(IN) :: ptimestep ! timestep (s)
+      REAL, INTENT(IN) :: fm(ngridmx,nlayermx+1) ! upward mass flux
+      REAL, INTENT(IN) :: entr(ngridmx,nlayermx) ! entrainment mass flux
+      REAL, INTENT(IN) :: detr(ngridmx,nlayermx) ! detrainment mass flux
+      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx) ! initial profil of q 
+      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx) ! mass of cells
+      INTEGER, INTENT(IN) :: ndt ! number of subtimesteps
+      INTEGER, INTENT(IN) :: zlmax ! index of maximum layer
 
 ! ============================ OUTPUTS ===========================
@@ -44,8 +48,9 @@
 ! =========== Init ==============================================
 
-      qa(:,:)=q_therm(:,:)
-      q(:,:)=q_therm(:,:)
+      qa(:,:)=q_therm(:,:) !q profile in the updraft
+      q(:,:)=q_therm(:,:) !mean q profile
 
 ! ====== Computing q ============================================
+! Based on equation 14 in appendix 4.2
 
       dq_therm(:,:)=0.
@@ -53,33 +58,32 @@
       invflux0(:,:)=ztimestep/masse0(:,:)      
 
-      do i=1,ndt
+      do i=1,ndt !subtimestep loop
 
-      do ig=1,ngrid
-         qa(ig,1)=q(ig,1)
-      enddo
+        do ig=1,ngrid
+           qa(ig,1)=q(ig,1)
+       enddo
 
-      do k=2,zlmax
-         do ig=1,ngridmx
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
-     &         1.e-5*masse0(ig,k)) then
-         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)
-            endif
-         enddo
-      enddo
+        do k=2,zlmax
+           do ig=1,ngridmx
+              if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
+     &        1.e-5*masse0(ig,k)) then
+                 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)
+              endif
+           enddo
+        enddo
 
-      do k=1,zlmax
-           q(:,k)=q(:,k)+         &
-     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
+        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
+     &    *invflux0(:,k)
+        enddo
 
       enddo !of do i=1,ndt
 
 ! ====== Derivative ==============================================
-
 
          do k=1,zlmax
Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 1020)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 1028)
@@ -1,2 +1,14 @@
+!=======================================================================
+! THERMCELL_MAIN_MARS
+! Author : A Colaitis
+! 
+! This routine is called by calltherm_interface and is inside a sub-timestep
+! loop. It computes thermals properties from parametrized entrainment and
+! detrainment rate as well as the source profile.
+! Mass flux are then computed and temperature and CO2 MMR are transported.
+!
+! This routine is based on the terrestrial version thermcell_main of the 
+! LMDZ model, written by C. Rio and F. Hourdin 
+!=======================================================================
 !
 !
@@ -4,9 +16,9 @@
      &                  ,pplay,pplev,pphi,zlev,zlay  &
      &                  ,pu,pv,pt,pq,pq2  &
-     &                  ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj  &
+     &                  ,pdtadj,pdqadj  &
      &                  ,fm,entr,detr,lmax,zmax  &
      &                  ,r_aspect &
      &                  ,zw2,fraca &
-     &                  ,zpopsk,ztla,heatFlux,heatFlux_down &
+     &                  ,zpopsk,heatFlux,heatFlux_down &
      &                  ,buoyancyOut, buoyancyEst)
 
@@ -14,59 +26,57 @@
 
 !=======================================================================
-! Mars version of thermcell_main. Author : A Colaitis
-!=======================================================================
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "tracer.h"
-#include "callkeys.h"
-
+
+#include "callkeys.h" !contains logical values determining which subroutines and which options are used. 
+                      ! includes logical "tracer", which is .true. if tracers are present and to be transported
+                      ! includes logical "lwrite", which is .true. if one wants more verbose outputs as GCM runs
+#include "dimensions.h" !contains global GCM grid dimension informations
+#include "dimphys.h" !similar to dimensions.h, and based on it; includes
+                     ! "ngridmx": number of horizontal grid points
+                     ! "nlayermx": number of atmospheric layers
+#include "comcstfi.h" !contains physical constant values such as
+                      ! "g" : gravitational acceleration (m.s-2)
+                      ! "r" : recuced gas constant (J.K-1.mol-1)
+#include "tracer.h" !contains tracer information such as
+                    ! "nqmx": number of tracers
+                    ! "noms()": tracer name
 ! ============== INPUTS ==============
 
-      REAL, INTENT(IN) :: ptimestep,r_aspect
-      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)
+      REAL, INTENT(IN) :: ptimestep !subtimestep (s)
+      REAL, INTENT(IN) :: r_aspect !aspect ratio (see paragraph 45 of paper and appendix S4)
+      REAL, INTENT(IN) :: pt(ngridmx,nlayermx) !temperature (K)
+      REAL, INTENT(IN) :: pu(ngridmx,nlayermx) !u component of the wind (ms-1)
+      REAL, INTENT(IN) :: pv(ngridmx,nlayermx) !v component of the wind (ms-1)
+      REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx) !tracer concentration (kg/kg)
+      REAL, INTENT(IN) :: pq2(ngridmx,nlayermx) ! Turbulent Kinetic Energy
+      REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) !Pressure at the middle of the layers (Pa)
+      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) !intermediate pressure levels (Pa)
+      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) !Geopotential at the middle of the layers (m2s-2)
+      REAL, INTENT(IN) :: zlay(ngridmx,nlayermx) ! altitude at the middle of the layers
+      REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
 
 ! ============== OUTPUTS ==============
 
-      REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
-      REAL :: pduadj(ngridmx,nlayermx)
-      REAL :: pdvadj(ngridmx,nlayermx)
-      REAL :: pdqadj(ngridmx,nlayermx,nqmx)
-!      REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
-      REAL :: pdq2adj(ngridmx,nlayermx)
-      REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
-
-! Diagnostics
+! TEMPERATURE
+      REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) !temperature change from thermals dT/dt (K/s)
+
+! DIAGNOSTICS
+      REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1) ! vertical velocity (m/s)
       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
+      REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
+
+! ============== LOCAL ================
+      REAL :: pdqadj(ngridmx,nlayermx,nqmx) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop)
 
 ! dummy variables when output not needed :
 
-!      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
 
-
 ! ============== LOCAL ================
 
       INTEGER ig,k,l,ll,iq
       INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
-      REAL linter(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)
@@ -85,5 +95,4 @@
 
       REAL wmax(ngridmx)
-      REAL wmax_sec(ngridmx)
       REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
 
@@ -115,5 +124,5 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL activecell(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega,adalim
+      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega
       INTEGER tic
 
@@ -145,24 +154,4 @@
       REAL f_old,ddd0,eee0,ddd,eee,zzz
       REAL fomass_max,alphamax
-
-! =========================================
-
-! ============= DTETA VARIABLES ===========
-
-! rien : on prends la divergence du flux turbulent
-
-! =========================================
-
-! ============= DV2 VARIABLES =============
-!               not used for now
-
-      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
 
 ! =========================================
@@ -183,18 +172,31 @@
 
 !-----------------------------------------------------------------------
-!   initialisation:
+!   initialization:
 !   ---------------
 
-      entr(:,:)=0.
-      detr(:,:)=0.
-      fm(:,:)=0.
-!      zu(:,:)=pu(:,:)
-!      zv(:,:)=pv(:,:)
-      zhc(:,:)=pt(:,:)/zpopsk(:,:)
+      entr(:,:)=0. ! entrainment mass flux
+      detr(:,:)=0. ! detrainment mass flux
+      fm(:,:)=0. ! upward mass flux
+      zhc(:,:)=pt(:,:)/zpopsk(:,:) ! potential temperature
       ndt=1
 
 ! **********************************************************************
-! Taking into account vertical molar mass gradients
+! In order to take into account the effect of vertical molar mass
+! gradient on convection, we define modified theta that depends
+! on the mass mixing ratio of Co2 in the cell.
+! See for details:
+!
+! Forget, F. and Millour, E. et al. "Non condensable gas enrichment and depletion
+! in the martian polar regions", third international workshop on the Mars Atmosphere:
+! Modeling and Observations, 1447, 9106. year: 2008
+!
+! This is especially important for modelling polar convection.
 ! **********************************************************************
+
+      
+!.......................................................................
+!  Special treatment for co2 at firstcall: compute coefficients and
+!  get index of co2 tracer
+!.......................................................................
 
       if(firstcall) then
@@ -235,5 +237,6 @@
        end if
 
-
+!------------------------------------------------------------------------
+! where are the different quantities defined ?
 !------------------------------------------------------------------------
 !                       --------------------
@@ -255,18 +258,7 @@
 
 !-----------------------------------------------------------------------
-!   Calcul des altitudes des couches
+!   Densities at layer and layer interface (see above), mass:
 !-----------------------------------------------------------------------
 
-!      do l=2,nlayermx
-!         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
-!      enddo
-!         zlev(:,1)=0.
-!         zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
-
-!         zlay(:,:)=pphi(:,:)/g
-!-----------------------------------------------------------------------
-!   Calcul des densites
-!-----------------------------------------------------------------------
-
       rho(:,:)=pplay(:,:)/(r*pt(:,:))
 
@@ -274,9 +266,8 @@
 
       do l=2,nlayermx
-!         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
           rhobarz(:,l)=pplev(:,l)/(r*0.5*(pt(:,l)+pt(:,l-1)))
       enddo
 
-!calcul de la masse
+! mass computation
       do l=1,nlayermx
          masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
@@ -284,4 +275,6 @@
 
 
+!-----------------------------------------------------------------
+!   Schematic representation of an updraft:
 !------------------------------------------------------------------
 !
@@ -330,11 +323,12 @@
 
 !=============================================================================
-!  Calculs initiaux ne faisant pas intervenir les changements de phase
+! Mars version: no phase change is considered, we use a "dry" definition
+! for the potential temperature.
 !=============================================================================
 
 !------------------------------------------------------------------
-!  1. alim_star est le profil vertical de l'alimentation a la base du
-!     panache thermique, calcule a partir de la flotabilite de l'air sec
-!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
+!  1. alim_star is the source layer vertical profile in the lowest layers
+! of the thermal plume. Computed from the air buoyancy
+!  2. lmin and lalim are the indices of begining and end of source profile
 !------------------------------------------------------------------
 !
@@ -343,12 +337,9 @@
 
 !-----------------------------------------------------------------------------
-!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
-!     panache sec conservatif (e=d=0) alimente selon alim_star 
-!     Il s'agit d'un calcul de type CAPE
-!     zmax_sec est utilise pour determiner la geometrie du thermique.
-!------------------------------------------------------------------------------
-!---------------------------------------------------------------------------------
-!calcul du melange et des variables dans le thermique
-!--------------------------------------------------------------------------------
+!  3. wmax and zmax are maximum vertical velocity and altitude of a 
+!     conservative plume (entrainment = detrainment = 0) using only
+!     the source layer. This is a CAPE computation used for determining 
+!     the closure mass flux. 
+!-----------------------------------------------------------------------------
 
 ! ===========================================================================
@@ -356,123 +347,56 @@
 ! ===========================================================================
 
-! Initialisations des variables reeles
-      ztva(:,:)=ztv(:,:)
-      ztva_est(:,:)=ztva(:,:)
-      ztla(:,:)=0.
-      zdz=0.
-      zbuoy(:,:)=0.
-      w_est(:,:)=0.
-      f_star(:,:)=0.
-      wa_moy(:,:)=0.
-      linter(:)=1.
+! Initialization
+      ztva(:,:)=ztv(:,:) ! temperature in the updraft = temperature of the env.
+      ztva_est(:,:)=ztva(:,:) ! estimated temp. in the updraft
+      ztla(:,:)=0. !intermediary variable
+      zdz=0. !layer thickness
+      zbuoy(:,:)=0. !buoyancy
+      w_est(:,:)=0. !estimated vertical velocity
+      f_star(:,:)=0. !non-dimensional upward mass flux f*
+      wa_moy(:,:)=0. !vertical velocity
 
 ! --------------------------------------------------------------------------
 ! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------
-! --------------  see thermiques.pro and getfit.py -------------------------
-
-!      a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6  ! svn baseline
-
-! Using broad downdraft selection
-!      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57
-!      ad = 0.0005114  ; bd = -0.662
-!      fdfu = -1.9
-
-! Using conditional sampling downdraft selection
-      a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421
-      ad = 0.00048088  ; bd = -0.6697
-!      fdfu = -1.3
-      fdfu=-0.8
-      a1inv=a1
-      b1inv=b1
-      omega=0.
-      adalim=0.
-
-! One good config for 34/35 levels
-!      a1inv=a1
-!      b1inv=b1
-!      be=1.1*be
-
-! Best configuration for 222 levels:
-
-!      omega=0.06
-!      b1=0.
-!      a1=1.
-!      a1inv=0.25*a1
-!      b1inv=0.0002
-!!
-!!
-!!      ae=0.9*ae
-
-! Best config for norad 222 levels:
-! and more specifically to C.large :
-
-!       omega=0.06
-!       omega=0.
-       omega=-0.03
-!       omega=0.
-       a1=1.
-!       b1=0.
-       b1=0.0001
-       a1inv=a1
-       be=1.1*be
-       ad = 0.0004
-!       ad=0.0003
-!       adalim=0.
-
-!       b1inv=0.00025
-       b1inv=b1
-
-!       b1inv = 0.0003
-
-!      omega=0.06
-! Trying stuff :
-
-!      ad=0.00035
-!      ae=0.95*ae
-!      b1=0.00055
-!      omega=0.04
-!
-!      ad = 0.0003
-!      ae=0.9*ae
-
-!      omega=0.04
-!!      b1=0.
-!      a1=1.
-!      a1inv=a1
-!      b1inv=0.0005689
-!!      be=1.1*be
-!!      ae=0.96*ae
-
-
-!       omega=0.06
-!       a1=1.
-!       b1=0.
-!       be=be
-!       a1inv=0.25*a1
-!       b1inv=0.0002    
-!       ad=1.1*ad
-!       ae=1.*ae
+
+! Detrainment
+      ad = 0.0004 !D_2 in paper, see paragraph 44
+      bd = -0.6697 !D_1 in paper, see paragraph 44
+
+! Entrainment
+      ae = 0.03683 !E_1 in paper, see paragraph 43
+      be = 0.631631 !E_2 in paper, see paragraph 43
+
+! Downdraft
+      fdfu=-0.8 !downdraft to updraft mass flux ratio, see paper paragraph 48
+      omega=-0.03 !see paper paragraph 48
+
+! Vertical velocity equation
+      a1=1. !a in paper, see paragraph 41
+      b1=0.0001 !b in paper, see paragraph 41
+
+! Inversion layer
+      a1inv=a1 !a1 in inversion layer
+      b1inv=b1 !b1 in inversion layer
 ! --------------------------------------------------------------------------
 ! --------------------------------------------------------------------------
 ! --------------------------------------------------------------------------
 
-! Initialisation des variables entieres
+! Some more initializations
       wmaxa(:)=0.
       lalim(:)=1
 
 !-------------------------------------------------------------------------
-! On ne considere comme actif que les colonnes dont les deux premieres
-! couches sont instables.
+! We consider as an activecell columns where the two first layers are
+! convectively unstable
+! When it is the case, we compute the source layer profile (alim_star)
+! see paper appendix 4.1 for details on the source layer
 !-------------------------------------------------------------------------
+
       activecell(:)=ztv(:,1)>ztv(:,2)
           do ig=1,ngridmx
             if (ztv(ig,1)>=(ztv(ig,2))) then
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
-!     &                       *log(1.+zlev(ig,2))
      &                       *sqrt(zlev(ig,2))
-!     &                       *sqrt(sqrt(zlev(ig,2)))
-!     &                       /sqrt(zlev(ig,2))
-!      &                       *zlev(ig,2)
-!      &                     *exp(-zlev(ig,2)/1000.)
                lalim(ig)=2
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
@@ -481,14 +405,9 @@
 
        do l=2,nlayermx-1
-!        do l=2,4
-         do ig=1,ngridmx
-           if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1).ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then
+         do ig=1,ngridmx
+           if (ztv(ig,l)>(ztv(ig,l+1)) .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.)  &
-!     &                       *log(1.+zlev(ig,l+1))
      &                       *sqrt(zlev(ig,l+1))
-!     &                       *sqrt(sqrt(zlev(ig,l+1)))
-!     &                       /sqrt(zlev(ig,l+1))
-!      &                       *zlev(ig,l+1)
-!      &                     *exp(-zlev(ig,l+1)/1000.)
                 lalim(ig)=l+1
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
@@ -505,45 +424,32 @@
 
       alim_star_tot(:)=1.
-!      if(alim_star(1,1) .ne. 0.) then
-!      print*, 'lalim star' ,lalim(1)
-!      endif
-
-!------------------------------------------------------------------------------
-! Calcul dans la premiere couche
-! On decide dans cette version que le thermique n'est actif que si la premiere
-! couche est instable.
-! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
-! dans une couche l>1
-!------------------------------------------------------------------------------
-
-      do ig=1,ngridmx
-! Le panache va prendre au debut les caracteristiques de l'air contenu
-! dans cette couche.
+
+! We compute the initial squared velocity (zw2) and non-dimensional upward mass flux
+! (f_star) in the first and second layer from the source profile.
+
+      do ig=1,ngridmx
           if (activecell(ig)) then
           ztla(ig,1)=ztv(ig,1)
-!cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1)
-! dans un panache conservatif
           f_star(ig,1)=0.
-          
           f_star(ig,2)=alim_star(ig,1)
-
           zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
       &                     *(zlev(ig,2)-zlev(ig,1))  &
-      &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
+      &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) !0.4=von Karman constant
           w_est(ig,2)=zw2(ig,2)
-
           endif
       enddo
 
-
 !==============================================================================
-!boucle de calcul de la vitesse verticale dans le thermique
+!==============================================================================
+!==============================================================================
+! LOOP ON VERTICAL LEVELS
 !==============================================================================
       do l=2,nlayermx-1
 !==============================================================================
-
-
-! On decide si le thermique est encore actif ou non
-! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
+!==============================================================================
+!==============================================================================
+
+
+! is the thermal plume still active ?
           do ig=1,ngridmx
              activecell(ig)=activecell(ig) &
@@ -553,25 +459,24 @@
 
 !---------------------------------------------------------------------------
-! calcul des proprietes thermodynamiques et de la vitesse de la couche l
-! sans tenir compte du detrainement et de l'entrainement dans cette
-! couche
-! C'est a dire qu'on suppose
-! ztla(l)=ztla(l-1)
-! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
-! avant) a l'alimentation pour avoir un calcul plus propre
+!
+! .I. INITIALIZATION
+!
+! Computations of the temperature and buoyancy properties in layer l, 
+! without accounting for entrainment and detrainment. We are therefore
+! assuming constant temperature in the updraft
+!
+! This computation yields an estimation of the buoyancy (zbuoy) and thereforce
+! an estimation of the velocity squared (w_est)
 !---------------------------------------------------------------------------
 
           do ig=1,ngridmx
              if(activecell(ig)) then
-!                if(l .lt. lalim(ig)) then
-!               ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
-!     &            alim_star(ig,l)*ztv(ig,l))  &
-!     &            /(f_star(ig,l)+alim_star(ig,l))
-!                else
                 ztva_est(ig,l)=ztla(ig,l-1)
-!                endif
 
                 zdz=zlev(ig,l+1)-zlev(ig,l)
                 zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+
+                ! Estimated vertical velocity squared 
+                ! (discretized version of equation 12 in paragraph 40 of paper)
 
                 if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
@@ -588,5 +493,5 @@
 
 !-------------------------------------------------
-!calcul des taux d'entrainement et de detrainement
+! Compute corresponding non-dimensional (ND) entrainment and detrainment rates
 !-------------------------------------------------
 
@@ -597,7 +502,9 @@
 
           if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
+
+         ! ND entrainment rate, see equation 16 of paper (paragraph 43)
+
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
         &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
-!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -607,29 +514,11 @@
              if(l .lt. lalim(ig)) then
 
-!                detr_star(ig,l)=0.
-                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
-            &  adalim
+                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)
-!                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
-!              &     0.0085*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.55)
-
-! last baseline from direct les
-!                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
-!               &     0.065*(2.5*(zbuoy(ig,l)/zw2m))**0.75
-
-! new param from continuity eq with a fit on dfdz
+         ! ND detrainment rate, see paragraph 44 of paper
+
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
             &  ad
-!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
-
-!               &     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)
-!                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
-!              &     ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002)
 
              endif
@@ -637,23 +526,9 @@
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
                 &     MAX(ad,bd*zbuoy(ig,l)/zw2m)
-!             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
-!             &  Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
-
-
-!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
-
-!              &  *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m)
-!              &  *5.*(-afact*zbuoy(ig,l)/zw2m)
-
-! last baseline from direct les
-!               &     0.065*(-2.5*(zbuoy(ig,l)/zw2m))**0.75
-
-! new param from continuity eq with a fit on dfdz
-
 
           endif
 
-! En dessous de lalim, on prend le max de alim_star et entr_star pour
-! alim_star et 0 sinon
+! If we are still in the source layer, we define the source layer entr. rate  (alim_star) as the
+! maximum between the source entrainment rate and the estimated entrainment rate.
 
        if (l.lt.lalim(ig)) then
@@ -662,5 +537,6 @@
        endif
 
-! Calcul du flux montant normalise
+! Compute the non-dimensional upward mass flux at layer l+1
+! using equation 11 of appendix 4.2 in paper
 
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
@@ -670,11 +546,31 @@
       enddo
 
-
-!----------------------------------------------------------------------------
-!calcul de la vitesse verticale en melangeant Tl et qt du thermique
-!---------------------------------------------------------------------------
-
+! -----------------------------------------------------------------------------------
+!
+! .II. CONVERGENCE LOOP
+!
+! We have estimated a vertical velocity profile and refined the source layer profile 
+! We now conduct iterations to compute:
+!
+! - the temperature inside the updraft from the estimated entrainment/source, detrainment,
+! and upward mass flux.
+! - the buoyancy from the new temperature inside the updraft
+! - the vertical velocity from the new buoyancy
+! - the entr., detr. and upward mass flux from the new buoyancy and vertical velocity
+!
+! This loop (tic) converges quickly. We have hardcoded 6 iterations from empirical observations.
+! Convergence occurs in 1 or 2 iterations in most cases.
+! -----------------------------------------------------------------------------------
+
+! -----------------------------------------------------------------------------------
+! -----------------------------------------------------------------------------------
       DO tic=0,5  ! internal convergence loop
+! -----------------------------------------------------------------------------------
+! -----------------------------------------------------------------------------------
+
+! Is the cell still active ?
       activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
+
+! If the cell is active, compute temperature inside updraft
       do ig=1,ngridmx
        if (activetmp(ig)) then
@@ -687,6 +583,8 @@
       enddo
 
+! Is the cell still active with respect to temperature variations ?
       activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
 
+! Compute new buoyancu and vertical velocity
       do ig=1,ngridmx
       if (activetmp(ig)) then
@@ -694,4 +592,5 @@
            zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
 
+          ! (discretized version of equation 12 in paragraph 40 of paper)
            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)-         &
@@ -704,4 +603,6 @@
 
 ! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
+         ! ND entrainment rate, see equation 16 of paper (paragraph 43)
+         ! ND detrainment rate, see paragraph 44 of paper
 
       do ig=1,ngridmx
@@ -713,5 +614,4 @@
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
         &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
-!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -721,12 +621,9 @@
              if(l .lt. lalim(ig)) then
 
-!                detr_star(ig,l)=0.
-                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
-            &  adalim
+                detr_star(ig,l)=0.
 
              else
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
             &  ad
-!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
 
              endif
@@ -734,5 +631,4 @@
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
                 &     MAX(ad,bd*zbuoy(ig,l)/zw2m)
-!             &  Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
 
           endif
@@ -742,6 +638,6 @@
           endif
 
-! En dessous de lalim, on prend le max de alim_star et entr_star pour
-! alim_star et 0 sinon
+! If we are still in the source layer, we define the source layer entr. rate  (alim_star) as the
+! maximum between the source entrainment rate and the estimated entrainment rate.
 
         if (l.lt.lalim(ig)) then
@@ -750,5 +646,6 @@
         endif
 
-! Calcul du flux montant normalise
+! Compute the non-dimensional upward mass flux at layer l+1
+! using equation 11 of appendix 4.2 in paper
 
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
@@ -757,9 +654,12 @@
       endif
       enddo
- 
-      ENDDO   ! of tic
+! -----------------------------------------------------------------------------------
+! -----------------------------------------------------------------------------------
+      ENDDO   ! of internal convergence loop
+! -----------------------------------------------------------------------------------
+! -----------------------------------------------------------------------------------
 
 !---------------------------------------------------------------------------
-!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
+! Miscellaneous computations for height 
 !---------------------------------------------------------------------------
 
@@ -767,13 +667,10 @@
             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
       IF (lwrite) THEN
-             print*,'On tombe sur le cas particulier de thermcell_plume'
+           print*,'thermcell_plume, particular case in velocity profile'
       ENDIF
                 zw2(ig,l+1)=0.
-                linter(ig)=l+1
             endif
 
         if (zw2(ig,l+1).lt.0.) then
-           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
-     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
            zw2(ig,l+1)=0.
         endif
@@ -781,5 +678,4 @@
 
         if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
-!   lmix est le niveau de la couche ou w (wa_moy) est maximum
             wmaxa(ig)=wa_moy(ig,l+1)
         endif
@@ -787,9 +683,16 @@
 
 !=========================================================================
-! FIN DE LA BOUCLE VERTICALE
-      enddo
 !=========================================================================
-
-!on recalcule alim_star_tot
+!=========================================================================
+! END OF THE LOOP ON VERTICAL LEVELS
+      enddo
+!=========================================================================
+!=========================================================================
+!=========================================================================
+
+! Recompute the source layer total entrainment alim_star_tot
+! as alim_star may have been modified in the above loop. Renormalization of
+! alim_star.
+
        do ig=1,ngridmx
           alim_star_tot(ig)=0.
@@ -817,11 +720,11 @@
 ! ===========================================================================
 
-! Attention, w2 est transforme en sa racine carree dans cette routine
+! WARNING, W2 (squared velocity) IS TRANSFORMED IN ITS SQUARE ROOT HERE
 
 !-------------------------------------------------------------------------------
-! Calcul des caracteristiques du thermique:zmax,wmax
+! Computations of the thermal height zmax and maximum vertical velocity wmax
 !-------------------------------------------------------------------------------
 
-!calcul de la hauteur max du thermique
+! Index of the thermal plume height
       do ig=1,ngridmx
          lmax(ig)=lalim(ig)
@@ -835,23 +738,13 @@
       enddo
 
-! On traite le cas particulier qu'il faudrait éviter ou le thermique
-! atteind le haut du modele ...
+! Particular case when the thermal reached the model top, which is not a good sign
       do ig=1,ngridmx
       if ( zw2(ig,nlayermx) > 1.e-10 ) then
-          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
+          print*,'WARNING !!!!! W2 non-zero in last layer'
           lmax(ig)=nlayermx
       endif
       enddo
 
-! pas de thermique si couche 1 stable
-!      do ig=1,ngridmx
-!         if (lmin(ig).gt.1) then
-!             lmax(ig)=1
-!             lmin(ig)=1
-!             lalim(ig)=1
-!         endif
-!      enddo
-!
-! Determination de zw2 max
+! Maximum vertical velocity zw2
       do ig=1,ngridmx
          wmax(ig)=0.
@@ -872,5 +765,8 @@
           enddo
       enddo
-!   Longueur caracteristique correspondant a la hauteur des thermiques.
+
+! Height of the thermal plume, defined as the following:
+! zmax=Integral[z*w(z)*dz]/Integral[w(z)*dz]
+!
       do  ig=1,ngridmx
          zmax(ig)=0.
@@ -892,6 +788,4 @@
        enddo
 
-! Attention, w2 est transforme en sa racine carree dans cette routine
-
 ! ===========================================================================
 ! ================= FIN HEIGHT ==============================================
@@ -904,5 +798,7 @@
       endif
 
-! Choix de la fonction d'alimentation utilisee pour la fermeture.
+! alim_star_clos is the source profile used for closure. It consists of the
+! modified source profile in the source layers, and the entrainment profile 
+! above it.
 
       alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
@@ -913,7 +809,7 @@
 
 !-------------------------------------------------------------------------------
-! Fermeture,determination de f
+! Closure, determination of the upward mass flux
 !-------------------------------------------------------------------------------
-! Appel avec la version seche
+! Init.
 
       alim_star2(:)=0.
@@ -921,5 +817,5 @@
       f(:)=0.
 
-! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
+! llmax is the index of the heighest thermal in the simulation domain
       llmax=1
       do ig=1,ngridmx
@@ -927,7 +823,6 @@
       enddo
 
-
-! Calcul des integrales sur la verticale de alim_star et de
-!   alim_star^2/(rho dz)
+! Integral of a**2/(rho* Delta z), see equation 13 of appendix 4.2 in paper
+
       do k=1,llmax-1
          do ig=1,ngridmx
@@ -940,12 +835,8 @@
       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
+! Closure mass flux, equation 13 of appendix 4.2 in paper
+
       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(ig)*alim_star_tot_clos(ig)/  &
       &     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
@@ -964,10 +855,14 @@
 
 !-------------------------------------------------------------------------------
-!deduction des flux
+! With the closure mass flux, we can compute the entrainment, detrainment and
+! upward mass flux from the non-dimensional ones.
 !-------------------------------------------------------------------------------
 
-      fomass_max=0.8
-      alphamax=0.5
-
+      fomass_max=0.8 !maximum mass fraction of a cell that can go upward in an
+                     ! updraft
+      alphamax=0.5 !maximum updraft coverage in a cell
+
+
+!    these variables allow to follow corrections made to the mass flux when lwrite=.true.
       ncorecfm1=0
       ncorecfm2=0
@@ -981,5 +876,5 @@
 
 !-------------------------------------------------------------------------
-! Multiplication par le flux de masse issu de la femreture
+! Multiply by the closure mass flux
 !-------------------------------------------------------------------------
 
@@ -988,4 +883,6 @@
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
+
+! Reconstruct the updraft mass flux everywhere
 
       do l=1,zlmax
@@ -1002,26 +899,22 @@
       enddo
 
-! Test provisoire : pour comprendre pourquoi on corrige plein de fois
-! le cas fm6, on commence par regarder une premiere fois avant les
-! autres corrections.
-
-!      do l=1,nlayermx
-!         do ig=1,ngridmx
-!            if (detr(ig,l).gt.fm(ig,l)) then
-!               ncorecfm8=ncorecfm8+1
-!            endif
-!         enddo
-!      enddo
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH Version en cours de test;
-! par rapport a thermcell_flux, on fait une grande boucle sur "l"
-! et on modifie le flux avec tous les contrï¿½les appliques d'affilee
-! pour la meme couche
-! Momentanement, on duplique le calcule du flux pour pouvoir comparer
-! les flux avant et apres modif
+!
+! Now we will reconstruct once again the upward
+! mass flux, but we will apply corrections
+! in some cases. We can compare to the
+! previously computed mass flux (above)
+!
+! This verification is done level by level
+!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      do l=1,zlmax
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      do l=1,zlmax !loop on the levels
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Upward mass flux at level l+1
 
          do ig=1,ngridmx
@@ -1038,5 +931,5 @@
 
 !-------------------------------------------------------------------------
-! Verification de la positivite des flux de masse
+! Upward mass flux should be positive
 !-------------------------------------------------------------------------
 
@@ -1061,42 +954,6 @@
          enddo
 
-!  Les "optimisations" du flux sont desactivecelles : moins de bidouilles
-!  je considere que celles ci ne sont pas justifiees ou trop delicates
-!  pour MARS, d'apres des observations LES. 
 !-------------------------------------------------------------------------
-!Test sur fraca croissant
-!-------------------------------------------------------------------------
-!      if (iflag_thermals_optflux==0) then
-!         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
-!!  zzz est le flux en l+1 a frac constant
-!             zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
-!     &                          /(rhobarz(ig,l)*zw2(ig,l))
-!             if (fm(ig,l+1).gt.zzz) then
-!                detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
-!                fm(ig,l+1)=zzz
-!                ncorecfm4=ncorecfm4+1
-!             endif
-!          endif
-!        enddo
-!      endif
-!
-!-------------------------------------------------------------------------
-!test sur flux de masse croissant
-!-------------------------------------------------------------------------
-!      if (iflag_thermals_optflux==0) then
-!         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)
-!              fm(ig,l+1)=fm(ig,l)
-!              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
-!               ncorecfm5=ncorecfm5+1
-!            endif
-!         enddo
-!      endif
-!
-!-------------------------------------------------------------------------
-!detr ne peut pas etre superieur a fm
+! Detrainment should be lower than upward mass flux
 !-------------------------------------------------------------------------
 
@@ -1107,10 +964,8 @@
                entr(ig,l)=fm(ig,l+1)
 
-! Dans le cas ou on est au dessus de la couche d'alimentation et que le
-! detrainement est plus fort que le flux de masse, on stope le thermique.
-!            endif
+! When detrainment is stronger than upward mass flux, and we are above the
+! thermal last level, the plume is stopped
 
             if(l.gt.lmax(ig)) then
-!            if(l.gt.lalim(ig)) then
                detr(ig,l)=0.
                fm(ig,l+1)=0.
@@ -1123,5 +978,5 @@
 
 !-------------------------------------------------------------------------
-!fm ne peut pas etre negatif
+! Check again for mass flux positivity
 !-------------------------------------------------------------------------
 
@@ -1135,19 +990,6 @@
 
 !-----------------------------------------------------------------------
-!la fraction couverte ne peut pas etre superieure a 1
+! Fractional coverage should be less than 1
 !-----------------------------------------------------------------------
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH Partie a revisiter.
-! Il semble qu'etaient codees ici deux optiques dans le cas
-! F/ (rho *w) > 1
-! soit limiter la hauteur du thermique en considerant que c'est
-! la derniere chouche, soit limiter F a rho w.
-! Dans le second cas, il faut en fait limiter a un peu moins
-! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
-! dans thermcell_main et qu'il semble de toutes facons deraisonable
-! d'avoir des fractions de 1..
-! Ci dessous, et dans l'etat actuel, le premier des  deux if est
-! sans doute inutile.
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
         do ig=1,ngridmx
@@ -1164,10 +1006,11 @@
         enddo
 
-! Fin de la grande boucle sur les niveaux verticaux
-      enddo
+      enddo ! on vertical levels
 
 !-----------------------------------------------------------------------
-! On fait en sorte que la quantite totale d'air entraine dans le
-! panache ne soit pas trop grande comparee a la masse de la maille
+!
+! We limit the total mass going from one level to the next, compared to the
+! initial total mass fo the cell
+!
 !-----------------------------------------------------------------------
 
@@ -1182,10 +1025,9 @@
                 entr(ig,l)=entr(ig,l)-eee
                 if ( ddd.gt.0.) then
-!   l'entrainement est trop fort mais l'exces peut etre compense par une
-!   diminution du detrainement)
+! The entrainment is too strong but we can compensate the excess by a detrainment decrease
                    detr(ig,l)=ddd
                 else
-!   l'entrainement est trop fort mais l'exces doit etre compense en partie
-!   par un entrainement plus fort dans la couche superieure
+! The entrainment is too strong and we compensate the excess by a stronger entrainment
+! in the layer above
                    if(l.eq.lmax(ig)) then
                       detr(ig,l)=fm(ig,l)+entr(ig,l)
@@ -1200,7 +1042,6 @@
          enddo
       enddo
-!
-!              ddd=detr(ig)-entre
-!on s assure que tout s annule bien en zmax
+
+! Check again that everything cancels at zmax
       do ig=1,ngridmx
          fm(ig,lmax(ig)+1)=0.
@@ -1210,5 +1051,6 @@
 
 !-----------------------------------------------------------------------
-! Impression du nombre de bidouilles qui ont ete necessaires
+! Summary of the number of modifications that were necessary (if lwrite=.true.
+! and only if there were a lot of them)
 !-----------------------------------------------------------------------
 
@@ -1237,19 +1079,13 @@
 
 !------------------------------------------------------------------
-!   calcul du transport vertical
+! vertical transport computation
 !------------------------------------------------------------------
 
 ! ------------------------------------------------------------------
-! Transport de teta dans l'updraft : (remplace thermcell_dq)
+! IN THE UPDRAFT
 ! ------------------------------------------------------------------
 
       zdthladj(:,:)=0.
-
-      if (1 .eq. 0) then
-!      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
-!     &     ,fm,entr,  &
-!     &    masse,ztv,zdthladj)
-      else
-
+! Based on equation 14 in appendix 4.2
 
       do ig=1,ngridmx
@@ -1270,8 +1106,6 @@
       enddo
 
-      endif
-
 ! ------------------------------------------------------------------
-! Prescription des proprietes du downdraft
+! DOWNDRAFT PARAMETERIZATION
 ! ------------------------------------------------------------------
 
@@ -1281,50 +1115,14 @@
          if (lmax(ig) .gt. 1) then
          do l=1,lmax(ig)
-!              if(zlay(ig,l) .le. 0.8*zmax(ig)) then
               if(zlay(ig,l) .le. zmax(ig)) then
+
+! see equation 18 of paragraph 48 in paper
               fm_down(ig,l) =fm(ig,l)* &
      &      max(fdfu,-4*max(0.,(zlay(ig,l)/zmax(ig)))-0.6)
               endif
 
-!             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
-!             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
-!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
-!             else
-!          ztvd(ig,l)=ztv(ig,l)
-!            endif
-
-!            if(zlay(ig,l) .le. 0.6*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
-!            elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
-!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
-!             else
-!          ztvd(ig,l)=ztv(ig,l)
-!            endif
-
-
-!            if(zlay(ig,l) .le. 0.8*zmax(ig)) then
-!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
-!            elseif(zlay(ig,l) .le. zmax(ig)) then
-!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
-!             else
-!          ztvd(ig,l)=ztv(ig,l)
-!            endif
-
-
-!             if (zbuoy(ig,l) .gt. 0.) then
-!               ztvd(ig,l)=ztva(ig,l)*0.9998
-!!               ztvd(ig,l)=ztv(ig,l)*0.997832
-!!             else
-!!               if(zlay(ig,l) .le. zmax(ig)) then            
-!!               ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
-!!               endif
-!             endif
-
             if(zlay(ig,l) .le. zmax(ig)) then            
+! see equation 19 of paragraph 49 in paper
           ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/400. + 0.997832))
-!          ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832))
              else
           ztvd(ig,l)=ztv(ig,l)
@@ -1336,5 +1134,5 @@
 
 ! ------------------------------------------------------------------
-! Transport de teta dans le downdraft : (remplace thermcell_dq)
+! TRANSPORT IN DOWNDRAFT
 ! ------------------------------------------------------------------
 
@@ -1344,4 +1142,5 @@
          if(lmax(ig) .gt. 1) then
 ! No downdraft in the very-near surface layer, we begin at k=3
+! Based on equation 14 in appendix 4.2
  
          do k=3,lmax(ig)
@@ -1361,5 +1160,5 @@
 
 !------------------------------------------------------------------
-! Calcul de la fraction de l'ascendance
+! Final fraction coverage with the clean upward mass flux, computed at interfaces
 !------------------------------------------------------------------
       fraca(:,:)=0.
@@ -1374,188 +1173,6 @@
       enddo
 
-
-
-! ===========================================================================
-! ============= DV2 =========================================================
-! ===========================================================================
-! ROUTINE OVERIDE : ne prends pas en compte le downdraft
-! de plus, le gradient de pression horizontal semble tout deregler... A VOIR
-
-      if (0 .eq. 1) then
-
 !------------------------------------------------------------------
-!  calcul du transport vertical du moment horizontal
-!------------------------------------------------------------------
-
-! Calcul du transport de V tenant compte d'echange par gradient
-! de pression horizontal avec l'environnement
-
-!   calcul du detrainement
-!---------------------------
-
-      nlarga0=0.
-
-      do k=1,nlayermx
-         do ig=1,ngridmx
-            detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
-         enddo
-      enddo
-
-!   calcul de la valeur dans les ascendances
-      do ig=1,ngridmx
-         zua(ig,1)=zu(ig,1)
-         zva(ig,1)=zv(ig,1)
-         ue(ig,1)=zu(ig,1)
-         ve(ig,1)=zv(ig,1)
-      enddo
-
-      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
-               gamma0(ig,k)=masse(ig,k)  &
-     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
-     &         *0.5/zmax(ig)  &
-     &         *1.
-            else
-               gamma0(ig,k)=0.
-            endif
-            if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1
-          enddo
-      enddo
-
-      gamma(:,:)=0.
-
-      do k=2,nlayermx
-
-         do ig=1,ngridmx
-
-            if (ltherm(ig,k)) then
-               dua(ig,k)=zua(ig,k-1)-zu(ig,k-1)
-               dva(ig,k)=zva(ig,k-1)-zv(ig,k-1)
-            else
-               zua(ig,k)=zu(ig,k)
-               zva(ig,k)=zv(ig,k)
-               ue(ig,k)=zu(ig,k)
-               ve(ig,k)=zv(ig,k)
-            endif
-         enddo
-
-
-! Debut des iterations
-!----------------------
-
-! AC WARNING : SALE !
-
-      do iter=1,5
-         do ig=1,ngridmx
-! Pour memoire : calcul prenant en compte la fraction reelle
-!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
-!              zf2=1./(1.-zf)
-! Calcul avec fraction infiniement petite
-               zf=0.
-               zf2=1.
-
-!  la première fois on multiplie le coefficient de freinage
-!  par le module du vent dans la couche en dessous.
-!  Mais pourquoi donc ???
-               if (ltherm(ig,k)) then
-!   On choisit une relaxation lineaire.
-!                 gamma(ig,k)=gamma0(ig,k)
-!   On choisit une relaxation quadratique.
-                gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
-                  zua(ig,k)=(fm(ig,k)*zua(ig,k-1)  &
-     &               +(zf2*entr(ig,k)+gamma(ig,k))*zu(ig,k))  &
-     &               /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2  &
-     &                 +gamma(ig,k))
-                  zva(ig,k)=(fm(ig,k)*zva(ig,k-1)  &
-     &               +(zf2*entr(ig,k)+gamma(ig,k))*zv(ig,k))  &
-     &               /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2  &
-     &                 +gamma(ig,k))
-
-!                  print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k)
-                  dua(ig,k)=zua(ig,k)-zu(ig,k)
-                  dva(ig,k)=zva(ig,k)-zv(ig,k)
-                  ue(ig,k)=(zu(ig,k)-zf*zua(ig,k))*zf2
-                  ve(ig,k)=(zv(ig,k)-zf*zva(ig,k))*zf2
-               endif
-         enddo
-! Fin des iterations
-!--------------------
-      enddo
-
-      enddo ! k=2,nlayermx
-
-! Calcul du flux vertical de moment dans l'environnement.
-!---------------------------------------------------------
-      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,ngridmx
-         wud(ig,1)=0.
-         wud(ig,nlayermx+1)=0.
-         wvd(ig,1)=0.
-         wvd(ig,nlayermx+1)=0.
-      enddo
-
-! calcul des tendances.
-!-----------------------
-      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)  &
-     &               -wud(ig,k)+wud(ig,k+1))  &
-     &               /masse(ig,k)
-            pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k)  &
-     &               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)  &
-     &               -wvd(ig,k)+wvd(ig,k+1))  &
-     &               /masse(ig,k)
-         enddo
-      enddo
-
-
-! Sorties eventuelles.
-!----------------------
-
-!      if (nlarga0>0) then
-!          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
-!      print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
-!          print*,'Il faudrait decortiquer ces points'
-!      endif
-
-! ===========================================================================
-! ============= FIN DV2 =====================================================
-! ===========================================================================
-
-      else
-!      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
-
-!------------------------------------------------------------------
-!  calcul du transport vertical de traceurs
+! Transport of C02 Tracer
 !------------------------------------------------------------------
 
@@ -1600,25 +1217,15 @@
          do ig=1,ngridmx
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
-!         pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l)
-         enddo
-      enddo
+         enddo
+      enddo
+
+! ===========================================================================
+! ============= FIN TRANSPORT ===============================================
+! ===========================================================================
+
 
 !------------------------------------------------------------------
-!  calcul du transport vertical de la tke
+!   Diagnostics for outputs
 !------------------------------------------------------------------
-
-!      modname='tke'
-!      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
-!      &      masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
-
-! ===========================================================================
-! ============= FIN TRANSPORT ===============================================
-! ===========================================================================
-
-
-!------------------------------------------------------------------
-!   Calculs de diagnostiques pour les sorties
-!------------------------------------------------------------------
-! DIAGNOSTIQUE
 ! We compute interface values for teta env and th. The last interface
 ! value does not matter, as the mass flux is 0 there.
