Index: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 506)
+++ trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90	(revision 508)
@@ -15,4 +15,5 @@
 #include "dimphys.h"
 #include "comcstfi.h"
+#include "tracer.h"
 
 !--------------------------------------------------------
@@ -61,5 +62,4 @@
       REAL ztla(ngridmx,nlayermx)
       REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
-      REAL dq_therm(ngridmx,nlayermx), dq_thermdown(ngridmx,nlayermx)
       REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
       REAL lmax_real(ngridmx)
@@ -100,6 +100,10 @@
       REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
 
-!---------------------------------------------------------
-
+!--------------------------------------------------------
+! Theta_m
+!--------------------------------------------------------
+
+      INTEGER ico2
+      SAVE ico2
 
 ! **********************************************************************
@@ -115,6 +119,4 @@
       q2_therm(:,:)=0.
       dq2_therm(:,:)=0.
-      dq_therm(:,:)=0.
-      dq_thermdown(:,:)=0.
       ztla(:,:)=0.
       pbl_dtke(:,:)=0.
@@ -165,4 +167,21 @@
 
 ! **********************************************************************
+! 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
+
+
+! **********************************************************************
 ! **********************************************************************
 ! **********************************************************************
@@ -204,6 +223,6 @@
          zheatFlux(:,:)=0.
          zheatFlux_down(:,:)=0.
-!         zbuoyancyOut(:,:)=0.
-!         zbuoyancyEst(:,:)=0.
+         zbuoyancyOut(:,:)=0.
+         zbuoyancyEst(:,:)=0.
          zzw2(:,:)=0.
          zmax(:)=0.
@@ -233,7 +252,7 @@
 !            d_v_the(:,:)=d_v_the(:,:)*fact
 !            dq2_the(:,:)=dq2_the(:,:)*fact
-!            if (nqmx .ne. 0) then
-!               d_q_the(:,:,:)=d_q_the(:,:,:)*fact
-!            endif
+            if (ico2 .ne. 0) then
+               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
+            endif
 
              zmaxth(:)=zmaxth(:)+zmax(:)*fact
@@ -250,8 +269,8 @@
             heatFlux_down(:,:)=heatFlux_down(:,:) &
      &       +zheatFlux_down(:,:)*fact
-!            buoyancyOut(:,:)=buoyancyOut(:,:) &
-!     &       +zbuoyancyOut(:,:)*fact
-!            buoyancyEst(:,:)=buoyancyEst(:,:) &
-!     &       +zbuoyancyEst(:,:)*fact
+            buoyancyOut(:,:)=buoyancyOut(:,:) &
+     &       +zbuoyancyOut(:,:)*fact
+            buoyancyEst(:,:)=buoyancyEst(:,:) &
+     &       +zbuoyancyEst(:,:)*fact
 
             zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
@@ -262,5 +281,7 @@
 !           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
 !           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
-!            d_q_ajs(:,:,:)=d_q_ajs(:,:,:)+d_q_the(:,:,:)
+            if (ico2 .ne. 0) then
+               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
+            endif
 !            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
 !  incrementation des variables meteo
@@ -269,5 +290,8 @@
 !            zu(:,:) = zu(:,:) + d_u_the(:,:)
 !            zv(:,:) = zv(:,:) + d_v_the(:,:)
-!            pq_therm(:,:,:) = pq_therm(:,:,:) + d_q_the(:,:,:)
+            if (ico2 .ne. 0) then
+             pq_therm(:,:,ico2) = &
+     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
+            endif
 !            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
 
@@ -303,8 +327,9 @@
       modname='tracer'
       DO iq=1,nqmx
+      if (iq .ne. ico2) then
       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
      &     ,fm_therm,entr_therm,detr_therm,  &
      &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
-
+      endif
       ENDDO
       endif
Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 506)
+++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90	(revision 508)
@@ -20,4 +20,6 @@
 #include "dimphys.h"
 #include "comcstfi.h"
+#include "tracer.h"
+#include "callkeys.h"
 
 ! ============== INPUTS ==============
@@ -110,5 +112,5 @@
       REAL wmaxa(ngridmx)
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
-      LOGICAL active(ngridmx),activetmp(ngridmx)
+      LOGICAL activecell(ngridmx),activetmp(ngridmx)
       REAL a1,b1,ae,be,ad,bd,fdfu
       INTEGER tic
@@ -163,4 +165,19 @@
 ! =========================================
 
+! ============== Theta_M Variables ========
+
+      INTEGER ico2
+      SAVE ico2
+      REAL m_co2, m_noco2, A , B
+      SAVE A, B
+      LOGICAL firstcall
+      save firstcall
+      data firstcall/.true./
+      REAL zhc(ngridmx,nlayermx)
+      REAL zdzfull(ngridmx,nlayermx)
+      REAL ratiom(ngridmx,nlayermx)
+
+! =========================================
+
 !-----------------------------------------------------------------------
 !   initialisation:
@@ -172,5 +189,47 @@
 !      zu(:,:)=pu(:,:)
 !      zv(:,:)=pv(:,:)
-      ztv(:,:)=pt(:,:)/zpopsk(:,:)
+      zhc(:,:)=pt(:,:)/zpopsk(:,:)
+
+! **********************************************************************
+! Taking into account vertical molar mass gradients
+! **********************************************************************
+
+      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
+                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
+                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
+!               Compute A and B coefficient use to compute
+!               mean molecular mass Mair defined by
+!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
+!               1/Mair = A*q(ico2) + B
+                A =(1/m_co2 - 1/m_noco2)
+                B=1/m_noco2
+             end if
+           enddo
+        endif
+
+      firstcall=.false.
+      endif !of if firstcall
+
+!.......................................................................
+!  Special treatment for co2
+!.......................................................................
+
+      if (ico2.ne.0) then
+!     Special case if one of the tracers is CO2 gas
+         DO l=1,nlayermx
+           DO ig=1,ngridmx
+            ztv(ig,l) = zhc(ig,l)*(A*pq(ig,l,ico2)+B)
+           ENDDO
+         ENDDO
+       else
+          ztv(:,:)=zhc(:,:)
+       end if
+
 
 !------------------------------------------------------------------------
@@ -332,5 +391,5 @@
 ! couches sont instables.
 !-------------------------------------------------------------------------
-      active(:)=ztv(:,1)>ztv(:,2)
+      activecell(:)=ztv(:,1)>ztv(:,2)
           do ig=1,ngridmx
             if (ztv(ig,1)>=(ztv(ig,2))) then
@@ -380,5 +439,5 @@
 ! Le panache va prendre au debut les caracteristiques de l'air contenu
 ! dans cette couche.
-          if (active(ig)) then
+          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)
@@ -407,5 +466,5 @@
 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
           do ig=1,ngridmx
-             active(ig)=active(ig) &
+             activecell(ig)=activecell(ig) &
       &                 .and. zw2(ig,l)>1.e-10 &
       &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
@@ -423,5 +482,5 @@
 
           do ig=1,ngridmx
-             if(active(ig)) then
+             if(activecell(ig)) then
 
 !                if(l .lt. lalim(ig)) then
@@ -453,5 +512,5 @@
 
       do ig=1,ngridmx
-        if (active(ig)) then
+        if (activecell(ig)) then
 
           zw2m=w_est(ig,l+1)
@@ -459,5 +518,6 @@
           if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
-        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+!        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+         &  MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -480,5 +540,6 @@
 ! new param from continuity eq with a fit on dfdz
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
-            &  ad
+!            &  ad
+             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
 
 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
@@ -492,5 +553,7 @@
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
-                &     bd*zbuoy(ig,l)/zw2m
+!                &     bd*zbuoy(ig,l)/zw2m
+             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+
 
 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
@@ -529,5 +592,5 @@
 
       DO tic=0,1  ! internal convergence loop
-      activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
+      activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
        if (activetmp(ig)) then
@@ -895,5 +958,5 @@
          enddo
 
-!  Les "optimisations" du flux sont desactivees : moins de bidouilles
+!  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. 
@@ -1332,4 +1395,35 @@
 
 !------------------------------------------------------------------
+!  calcul du transport vertical de traceurs
+!------------------------------------------------------------------
+
+! We only transport co2 tracer because it is coupled to the scheme through theta_m
+! The rest is transported outside the sub-timestep loop
+
+      if (ico2.ne.0) then
+!      if (nqmx .ne. 0.) then
+      do l=1,nlayermx
+         zdzfull(:,l)=zlev(:,l+1)-zlev(:,l)
+      enddo
+
+      modname='tracer'
+      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
+     &     ,fm,entr,detr,  &
+     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),modname,zdzfull)
+!      endif
+
+! Compute the ratio between theta and theta_m
+      
+       do l=1,nlayermx
+          do ig=1,ngridmx
+             ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
+          enddo
+       enddo
+       else
+          ratiom(:,:)=1.
+       endif
+
+
+!------------------------------------------------------------------
 !  incrementation dt
 !------------------------------------------------------------------
@@ -1337,21 +1431,7 @@
       do l=1,nlayermx
          do ig=1,ngridmx
-           pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
-!           pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)
-         enddo
-      enddo
-
-!------------------------------------------------------------------
-!  calcul du transport vertical de traceurs
-!------------------------------------------------------------------
-
-!      if (nqmx .ne. 0.) then
-!      modname='tracer'
-!      DO iq=1,nqmx
-!          call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
-!          &      masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax)
-!
-!      ENDDO
-!      endif
+         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
+         enddo
+      enddo
 
 !------------------------------------------------------------------
@@ -1378,7 +1458,7 @@
       do l=1,nlayermx-1
        do ig=1,ngridmx
-        teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))
-        teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))
-        teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))
+        teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))*ratiom(ig,l)
+        teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))*ratiom(ig,l)
+        teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))*ratiom(ig,l)
        enddo
       enddo
@@ -1391,6 +1471,6 @@
        do ig=1,ngridmx
         heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
-!        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-!        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
         heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
        enddo
