Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 472)
+++ trunk/LMDZ.MARS/README	(revision 473)
@@ -1339,2 +1339,12 @@
    order of tracers was changed).
 
+== 14/12/2011 == AC
+>> Improved computation of mixing in vdifc for temperatures less than the saturation
+   temperature.
+   The scheme works by computing mixing of T and then liberating latent heat when T<Ts.
+   Taking into account the mass of co2 ice condensated, a new mixing profile for T is computed.
+   When the computation has converged, the tendancy for mixing of T is taken with respect to
+   the initial T profile for which co2 ice has been condensated.
+>> This scheme is based on initial work done by Melanie Vangvichith and has been adapted/tested/bug-fixed
+   and improved to Mars.
+
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 472)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 473)
@@ -72,4 +72,6 @@
 c   ------
 
+      REAL ust(ngridmx),tst(ngridmx)
+ 
       INTEGER ilev,ig,ilay,nlev
 
@@ -112,4 +114,28 @@
       REAL kmixmin
 
+
+c    Mass-variation scheme :
+c    ~~~~~~~
+
+      INTEGER j,l
+      REAL zcondicea(ngrid,nlayermx)
+      REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1)
+      REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx)
+      REAL pdtc(ngrid,nlayermx)
+      REAL zhs(ngridmx,nlayermx)
+      REAL cpice,ccond
+      SAVE ccond,cpice
+      DATA cpice /1000./
+
+c     Theta_m formulation for mass-variation scheme :
+c     ~~~~~~~
+
+      INTEGER ico2,llnt(ngridmx)
+      SAVE ico2
+      REAL m_co2, m_noco2, A , B
+      SAVE A, B, m_co2, m_noco2
+      REAL vmr_co2(ngridmx,nlayermx)
+      REAL qco2,mmean
+
 c    ** un petit test de coherence
 c       --------------------------
@@ -126,6 +152,27 @@
          bcond=1./tcond1mb
          acond=r/latcond
+         ccond=cpp/(g*latcond)
          PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
-         PRINT*,'          acond,bcond',acond,bcond
+         PRINT*,'          acond,bcond,ccond',acond,bcond,ccond
+
+
+         ico2=0
+
+         if (tracer) then
+c          Prepare Special treatment if one of the tracer 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)
+c               Compute A and B coefficient use to compute
+c               mean molecular mass Mair defined by
+c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
+c               1/Mair = A*q(ico2) + B
+                A =(1/m_co2 - 1/m_noco2)
+                B=1/m_noco2
+             endif
+           enddo
+         end if
 
         firstcall=.false.
@@ -149,4 +196,6 @@
          DO ig=1,ngrid
             za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
+! Mass variation scheme:
+            betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay))
          ENDDO
       ENDDO
@@ -163,5 +212,5 @@
       ENDDO
       DO ig=1,ngrid
-		 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
+         zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
       ENDDO
 
@@ -185,17 +234,51 @@
       ENDIF
 
+c     -----------------------------------
 c     Potential Condensation temperature:
 c     -----------------------------------
 
-c     if (callcond) then 
-c       DO ilev=1,nlay
-c         DO ig=1,ngrid
-c           zhcond(ig,ilev) =
-c    &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
-c         END DO
-c       END DO
-c     else
-        call zerophys(ngrid*nlay,zhcond)
-c     end if
+c     Compute CO2 Volume mixing ratio
+c     -------------------------------
+      if (callcond.and.(ico2.ne.0)) then
+         DO ilev=1,nlay
+            DO ig=1,ngrid
+              qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep
+c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
+              mmean=1/(A*qco2 +B)
+              vmr_co2(ig,ilev) = qco2*mmean/m_co2
+            ENDDO
+         ENDDO
+      else
+         DO ilev=1,nlay
+            DO ig=1,ngrid
+              vmr_co2(ig,ilev)=0.95
+            ENDDO
+         ENDDO
+      end if
+
+c     forecast of atmospheric temperature zt and frost temperature ztcond
+c     --------------------------------------------------------------------
+
+      DO ilev=1,nlay
+         DO ig=1,ngrid
+            ztcond(ig,ilev)=
+     &         1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev)))
+            if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
+         ENDDO
+      ENDDO
+
+      ztcond(:,nlay+1)=ztcond(:,nlay)
+
+      if (callcond) then 
+        DO ilev=1,nlay
+          DO ig=1,ngrid
+!            zhcond(ig,ilev) =
+!     &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
+              zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
+          END DO
+        END DO
+      else
+         call zerophys(ngrid*nlay,zhcond)
+      end if
 
 
@@ -209,5 +292,5 @@
             zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
             zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
-            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
+!            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
          ENDDO
       ENDDO
@@ -239,7 +322,13 @@
         zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
         zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
+!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)
+!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
+!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
+!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
         ELSE
         zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance 
         zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
+!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
+!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
         ENDIF
 
@@ -254,11 +343,7 @@
 !     &                         2,zcdh_true)
 !        call WRITEDIAGFI(ngridmx,'u*',
-!     &              'friction velocity','m/s',
-!     &             2,sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2))
-!        call WRITEDIAGFI(ngridmx,'T*',
-!     &              'friction temperature','K',
-!     &   2,zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)*
-!     &   -(ptsrf(:)-ph(:,1))
-!     & /(sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)))
+!     &              'friction velocity','m/s',2,ust)
+!       call WRITEDIAGFI(ngridmx,'T*',
+!     &              'friction temperature','K',2,tst)
 !        call WRITEDIAGFI(ngridmx,'rm-1',
 !     &              'aerodyn momentum conductance','m/s',
@@ -405,22 +490,13 @@
 c       -------------
 
+c Mass variation scheme:
       CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
       CALL multipl(ngrid,zcdh,zb0,zb)
 
-      DO ig=1,ngrid
-         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
-         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-      ENDDO
-
-      DO ilay=nlay-1,1,-1
-         DO ig=1,ngrid
-            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
-     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-         ENDDO
-      ENDDO
+c on initialise dm c
+      
+      pdtc(:,:)=0.
+      zt(:,:)=0.
+      dmice(:,:)=0.
 
 c    ** calcul de (d Planck / dT) a la temperature d'interface
@@ -431,4 +507,57 @@
          zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
       ENDDO
+
+
+! calcul de zc et zd pour la couche top en prenant en compte le terme
+! de variation de masse (on fait une boucle pour que ça converge)
+
+! Identification des points de grilles qui ont besoin de la correction
+
+      llnt(:)=1
+
+      DO ig=1,ngrid
+         DO l=1,nlay
+            if(zh(ig,l) .lt. zhcond(ig,l)) then
+               llnt(ig)=300  
+! 200 and 100 do not go beyond month 9 with normal dissipation
+               goto 5
+            endif
+         ENDDO
+5     continue
+      ENDDO
+
+      DO ig=1,ngrid
+
+! Initialization of z1 and zd, which do not depend on dmice
+
+      z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+      zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+
+      DO ilay=nlay-1,1,-1
+          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $        zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+      ENDDO
+
+! Convergence loop
+
+      DO j=1,llnt(ig)
+
+            z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+            zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
+     &      -betam(ig,nlay)*dmice(ig,nlay)
+            zc(ig,nlay)=zc(ig,nlay)*z1(ig)
+!            zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+
+! calcul de zc et zd pour les couches du haut vers le bas
+
+             DO ilay=nlay-1,1,-1
+               z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+               zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
+     $            zb(ig,ilay+1)*zc(ig,ilay+1)-
+     $            betam(ig,ilay)*dmice(ig,ilay))*z1(ig)
+!               zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+            ENDDO
 
 c    ** calcul de la temperature_d'interface et de sa tendance.
@@ -444,30 +573,43 @@
 c       ----------------------------------------------------
 
-      DO ig=1,ngrid
          z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
      s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
          z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
          ztsrf2(ig)=z1(ig)/z2(ig)
-         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
-
-c        Modif speciale CO2 condensation:
-c        tconds = 1./(bcond-acond*log(.0095*pplev(ig,1)))
-c        if ((callcond).and.
-c    &      ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
-c           zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
-c        else
-            zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
-c        end if
-      ENDDO
+!         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep  !incremented outside loop
+
+            zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
 
 c    ** et a partir de la temperature au sol on remonte 
 c       -----------------------------------------------
 
-      DO ilay=2,nlay
-         DO ig=1,ngrid
-            hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond
-            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
-         ENDDO
-      ENDDO
+            DO ilay=2,nlay
+               zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
+            ENDDO
+            DO ilay=1,nlay
+               zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
+            ENDDO
+
+c      Condensation/sublimation in the atmosphere
+c      ------------------------------------------
+c      (computation of zcondicea and dmice)
+
+      zcondicea(ig,:)=0.
+      pdtc(ig,:)=0.
+
+      DO l=nlay , 1, -1
+           IF(zt(ig,l).LT.ztcond(ig,l)) THEN
+               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
+               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
+     &                        *ccond*pdtc(ig,l)
+              dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
+            END IF
+      ENDDO
+
+         ENDDO    !of Do j=1,XXX 
+
+      ENDDO   !of Do ig=1,nlay
+
+      pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
 
 
@@ -570,5 +712,5 @@
             ! special case for water ice tracer: do not include
             ! h2o ice tracer from surface (which is set when handling
-            ! h2o vapour case (see further down).
+           ! h2o vapour case (see further down).
             DO ig=1,ngrid
                 z1(ig)=1./(za(ig,1)+zb(ig,1)+
@@ -641,11 +783,10 @@
             pdvdif(ig,ilev)=(    zv(ig,ilev)-
      $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
-            hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
-     $           zhcond(ig,ilev))        ! modif co2cond
-            pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
-         ENDDO
-      ENDDO
-
-    
+            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 
+     $  + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
+            pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
+         ENDDO
+      ENDDO
+
       if (tracer) then 
         DO iq = 1, nq
@@ -670,18 +811,9 @@
          DO ilev=1,nlay
             WRITE(*,'(4f15.7)')
-     s      ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
+     s      ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
      s      pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
 
          ENDDO
       ENDIF
-
-      if (calltherm .and. outptherm) then
-      if (ngrid .eq. 1) then
-        call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
-     &                       'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
-        call WRITEDIAGFI(ngrid,'zkh','zkh',
-     &                       'SI',1,zkh)
-      endif
-      endif 
 
       RETURN
