Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 682)
+++ 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.
