Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2009)
+++ trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2009)
@@ -0,0 +1,766 @@
+      MODULE co2condens_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep,
+     $                  pcapcal,pplay,pplev,ptsrf,pt,
+     $                  pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
+     $                  piceco2,psolaralb,pemisurf,
+     $                  pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
+     $                  fluxsurf_sw,zls,
+     $                  zdqssed_co2,pcondicea_co2microp,
+     $                  zdtcloudco2)
+                                                   
+       use tracer_mod, only: noms
+       use surfdat_h, only: emissiv, phisfi
+       use geometry_mod, only: latitude ! grid point latitudes (rad)
+       use planete_h, only: obliquit
+       use comcstfi_h, only: cpp, g, r, pi 
+#ifndef MESOSCALE
+       USE vertical_layers_mod, ONLY: bp
+#endif
+       IMPLICIT NONE
+c=======================================================================
+c   subject:
+c   --------
+c   Condensation/sublimation of CO2 ice on the ground and in the
+c   atmosphere 
+c   (Scheme described in Forget et al., Icarus, 1998)
+c
+c   author:   Francois Forget     1994-1996 ; updated 1996 -- 2018
+c   ------
+c             adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
+c
+c=======================================================================
+c
+c    0.  Declarations :
+c    ------------------
+c
+      include "callkeys.h"
+
+c-----------------------------------------------------------------------
+c    Arguments :
+c    ---------
+      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
+      INTEGER,INTENT(IN) :: nq  ! number of tracers
+
+      REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
+      REAL,INTENT(IN) :: pcapcal(ngrid)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
+      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
+      REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
+                                           ! previous physical processes (K/s)
+      REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
+                                           ! from previous physical processes
+      REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
+                                           ! from previous physical processes
+      REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
+                                       ! previous physical processes (K/s)
+      REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
+      REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
+      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
+                                              ! previous physical processes
+
+      REAL,INTENT(IN) :: zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
+      REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
+      REAL,INTENT(IN) :: zdtcloudco2(ngrid,nlayer) ! tendency on temperature due to latent heat                 
+
+      REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
+      REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
+      REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
+      
+      ! tendencies due to CO2 condensation/sublimation:
+      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
+      REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
+      REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
+      REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
+      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
+      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
+      
+      ! added to calculate flux dependent albedo:
+      REAL,intent(in) :: fluxsurf_sw(ngrid,2)
+      real,intent(in) :: zls ! solar longitude (rad)
+
+c
+c    Local variables :
+c    -----------------
+
+      INTEGER i,j
+      INTEGER l,ig,iq,icap
+      REAL zt(ngrid,nlayer)
+      REAL zcpi
+      REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
+      REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
+      REAL zdiceco2(ngrid)
+      REAL zcondicea(ngrid,nlayer) ! condensation rate in layer  l  (kg/m2/s)
+      REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s)
+      REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
+      REAL zfallheat
+      REAL zmflux(nlayer+1)
+      REAL zu(nlayer),zv(nlayer)
+      REAL zq(nlayer,nq),zq1(nlayer)
+      REAL ztsrf(ngrid)
+      REAL ztc(nlayer), ztm(nlayer+1) 
+      REAL zum(nlayer+1) , zvm(nlayer+1)
+      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
+      REAL masse(nlayer),w(nlayer+1)
+      REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
+      LOGICAL condsub(ngrid)
+
+      real :: emisref(ngrid)
+
+c variable speciale diagnostique
+      real tconda1(ngrid,nlayer)
+      real tconda2(ngrid,nlayer)
+c     REAL zdiceco2a(ngrid) ! for diagnostic only
+      real zdtsig (ngrid,nlayer)
+      real zdt (ngrid,nlayer)
+      real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
+! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
+! then condensation temperature is computed using partial pressure of CO2
+      logical,parameter :: improved_ztcond=.true.
+
+c   local saved variables
+      integer,save :: ico2 ! index of CO2 tracer
+      real,save :: qco2min,qco2,mmean
+      real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
+      real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
+      real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
+      REAL,SAVE :: acond,bcond,ccond
+      real,save :: m_co2, m_noco2, A , B
+
+      LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
+
+c D.BARDET: to debug 
+      real ztc3D(ngrid,nlayer)
+      REAL ztm3D(ngrid,nlayer)
+      REAL zmflux3D(ngrid,nlayer)
+c----------------------------------------------------------------------
+
+c   Initialisation
+c   --------------
+c
+      ! AS: firstcall OK absolute
+      IF (firstcall) THEN
+         
+         bcond=1./tcond1mb
+         ccond=cpp/(g*latcond)
+         acond=r/latcond
+
+         firstcall=.false.
+         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
+         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
+         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,nq
+             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
+c          minimum CO2 mix. ratio below which mixing occur with layer above: 
+           qco2min =0.75  
+         end if
+      ENDIF ! of IF (firstcall)
+      zcpi=1./cpp
+
+c
+c======================================================================
+c    Calcul of CO2 condensation sublimation 
+c    ============================================================
+c  
+c    Used variable :
+c       piceco2(ngrid)   :  amount of co2 ice on the ground (kg/m2)
+c       zcondicea(ngrid,l):  condensation rate in layer  l  (kg/m2/s)
+c       zcondices(ngrid):  condensation rate on the ground (kg/m2/s)
+c       zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
+c                           
+c       pdtc(ngrid,nlayer)  : dT/dt due to cond/sub
+c
+c
+c     Tendencies set to 0 (except pdtc)
+c     -------------------------------------
+      DO l=1,nlayer
+         DO ig=1,ngrid
+           zcondicea(ig,l) = 0.
+           zfallice(ig,l) = 0.
+           pduc(ig,l)  = 0
+           pdvc(ig,l)  = 0
+         END DO
+      END DO
+         
+      DO iq=1,nq
+      DO l=1,nlayer
+         DO ig=1,ngrid
+           pdqc(ig,l,iq)  = 0
+         END DO
+      END DO
+      END DO
+
+      DO ig=1,ngrid
+         zfallice(ig,nlayer+1) = 0.
+         zcondices(ig) = 0.
+         pdtsrfc(ig) = 0.
+         pdpsrf(ig) = 0.
+         condsub(ig) = .false.
+         zdiceco2(ig) = 0.
+      ENDDO
+      zfallheat=0
+
+c     *************************
+c     ATMOSPHERIC CONDENSATION
+c     *************************
+
+c     Compute CO2 Volume mixing ratio
+c     -------------------------------
+      if (improved_ztcond.and.(ico2.ne.0)) then
+         DO l=1,nlayer
+            DO ig=1,ngrid
+              qco2=pq(ig,l,ico2)+pdq(ig,l,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,l) = qco2*mmean/m_co2 
+            ENDDO
+         ENDDO
+      else
+         DO l=1,nlayer
+            DO ig=1,ngrid
+              vmr_co2(ig,l)=0.95
+            ENDDO
+         ENDDO
+      end if
+
+      IF (.NOT. co2clouds) then
+c     forecast of atmospheric temperature zt and frost temperature ztcond
+c     --------------------------------------------------------------------
+
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
+!            ztcond(ig,l)=1./(bcond-acond*log(.0095*pplay(ig,l)))
+            if (pplay(ig,l).ge.1e-4) then
+              ztcond(ig,l)=
+     &         1./(bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l)))
+            else
+              ztcond(ig,l)=0.0 !mars Monica
+            endif
+         ENDDO
+      ENDDO
+
+      ztcond(:,nlayer+1)=ztcond(:,nlayer)
+ 
+c      Condensation/sublimation in the atmosphere
+c      ------------------------------------------
+c      (calcul of zcondicea , zfallice and pdtc)
+c
+      DO l=nlayer , 1, -1
+         DO ig=1,ngrid
+           pdtc(ig,l)=0.
+           IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
+               condsub(ig)=.true.
+               IF (zfallice(ig,l+1).gt.0) then  
+                 zfallheat=zfallice(ig,l+1)*
+     &           (pphi(ig,l+1)-pphi(ig,l) +
+     &          cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
+               ELSE
+                   zfallheat=0.
+               ENDIF
+               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
+               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
+     &                        *ccond*pdtc(ig,l)- zfallheat
+c              Case when the ice from above sublimes entirely
+c              """""""""""""""""""""""""""""""""""""""""""""""
+               IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
+                  pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
+     &              (ccond*(pplev(ig,l)-pplev(ig,l+1)))
+                  zcondicea(ig,l)= -zfallice(ig,l+1)
+               END IF
+
+               zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
+            END IF
+         ENDDO
+      ENDDO
+      
+      ELSE
+         
+        DO ig=1,ngrid 
+         zfallice(ig,1) = zdqssed_co2(ig)
+        ENDDO
+        DO l=nlayer , 1, -1
+            DO ig=1,ngrid 
+         zcondicea(ig,l) = pcondicea_co2microp(ig,l)* 
+     &                        (pplev(ig,l) - pplev(ig,l+1))/g
+            ENDDO
+        ENDDO
+        DO l=nlayer, 1, -1
+         DO ig=1, ngrid
+            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
+            pdtc(ig,l)=0.
+         ENDDO
+      ENDDO 
+      ENDIF ! if not co2clouds
+
+      call WRITEdiagfi(ngrid,"pdtc_atm",
+     &         "temperature tendency due to CO2 condensation",
+     &         " ",3,pdtc) 
+     
+       call WRITEdiagfi(ngrid,"zcondicea",
+     &         "",
+     &         " ",3,zcondicea)         
+      
+       call WRITEdiagfi(ngrid,"zfallice",
+     &         "",
+     &         " ",2,zfallice(ngrid,1))
+
+c     *************************
+c       SURFACE  CONDENSATION
+c     *************************
+
+c     forecast of ground temperature ztsrf and frost temperature ztcondsol
+c     --------------------------------------------------------------------
+      DO ig=1,ngrid
+         ztcondsol(ig)=
+     &          1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1)))
+         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
+      ENDDO
+
+c
+c      Condensation/sublimation on the ground 
+c      --------------------------------------
+c      (compute zcondices and pdtsrfc)
+c
+      DO ig=1,ngrid
+         IF(latitude(ig).lt.0) THEN
+           ! Southern hemisphere
+            icap=2
+         ELSE
+           ! Northern hemisphere
+            icap=1
+         ENDIF
+       
+c
+c        Loop on where we have  condensation/ sublimation 
+         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.   ! ground cond 
+     $       (zfallice(ig,1).NE.0.) .OR.           ! falling snow
+     $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
+     $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
+            condsub(ig) = .true. 
+
+            IF (zfallice(ig,1).gt.0) then  
+                 zfallheat=zfallice(ig,1)*
+     &           (pphi(ig,1)- phisfi(ig) +
+     &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
+            ELSE
+                 zfallheat=0.
+            ENDIF
+
+c           condensation or partial sublimation of CO2 ice
+c           """""""""""""""""""""""""""""""""""""""""""""""
+            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 
+     &      /(latcond*ptimestep)         - zfallheat
+            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
+       
+c           If the entire CO_2 ice layer sublimes
+c           """"""""""""""""""""""""""""""""""""""""""""""""""""
+c           (including what has just condensed in the atmosphere)
+
+            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
+     &          -zcondices(ig))THEN
+               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
+               pdtsrfc(ig)=(latcond/pcapcal(ig))*
+     &         (zcondices(ig)+zfallheat)
+            END IF
+
+c           Changing CO2 ice amount and pressure :
+c           """"""""""""""""""""""""""""""""""""
+
+            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
+            piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
+            pdpsrf(ig) = -zdiceco2(ig)*g
+
+            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
+               PRINT*,'STOP in condens'
+               PRINT*,'condensing more than total mass'
+               PRINT*,'Grid point ',ig
+               PRINT*,'Ps = ',pplev(ig,1)
+               PRINT*,'d Ps = ',pdpsrf(ig)
+               STOP
+            ENDIF
+         END IF ! if there is condensation/sublimmation
+      ENDDO ! of DO ig=1,ngrid
+
+c  ********************************************************************
+c  Surface albedo and emissivity of the surface below the snow (emisref)
+c  ********************************************************************
+
+!     Check that amont of CO2 ice is not problematic
+      DO ig=1,ngrid
+           if(.not.piceco2(ig).ge.0.) THEN
+              if(piceco2(ig).le.-5.e-8) print*,
+     $         'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
+               piceco2(ig)=0.
+           endif
+      ENDDO
+      
+!     Set albedo and emissivity of the surface
+!     ----------------------------------------
+      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
+
+! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
+      DO ig=1,ngrid
+        if (piceco2(ig).eq.0) then
+          pemisurf(ig)=emissiv
+        endif
+      ENDDO
+
+!         firstcall2=.false.
+c ***************************************************************
+c  Correction to account for redistribution between sigma or hybrid 
+c  layers when changing surface pressure (and warming/cooling
+c  of the CO2 currently changing phase).
+c *************************************************************
+
+      DO ig=1,ngrid
+        if (condsub(ig)) then
+           do l=1,nlayer
+             ztc(l)  =zt(ig,l)   +pdtc(ig,l)  *ptimestep
+             zu(l)   =pu(ig,l)   +pdu( ig,l)  *ptimestep
+             zv(l)   =pv(ig,l)   +pdv( ig,l)  *ptimestep
+            do iq=1,nq
+             zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
+            enddo
+           end do
+
+c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
+c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+            zmflux(1) = -zcondices(ig)
+            DO l=1,nlayer
+             zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
+#ifndef MESOSCALE
+     &        + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
+c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
+          if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
+#else
+          if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
+#endif
+            END DO
+
+c Mass of each layer
+c ------------------ 
+            DO l=1,nlayer
+              masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
+            END DO
+
+
+c  Corresponding fluxes for T,U,V,Q
+c  """"""""""""""""""""""""""""""""
+
+c           averaging operator for TRANSPORT  
+c           """"""""""""""""""""""""""""""""
+c           Value transfert at the surface interface when condensation
+c           sublimation:
+            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
+            zum(1) = 0  
+            zvm(1) = 0  
+            do iq=1,nq
+              zqm(1,iq)=0. ! most tracer do not condense !
+            enddo
+c           Special case if one of the tracer is CO2 gas
+            if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
+
+c           Van Leer scheme:
+            DO l=1,nlayer+1
+                w(l)=-zmflux(l)*ptimestep
+            END DO
+            call vl1d(nlayer,ztc,2.,masse,w,ztm) 
+            call vl1d(nlayer,zu ,2.,masse,w,zum) 
+            call vl1d(nlayer,zv ,2.,masse,w,zvm) 
+            do iq=1,nq
+             do l=1,nlayer
+              zq1(l)=zq(l,iq)
+             enddo
+             zqm1(1)=zqm(1,iq)
+             call vl1d(nlayer,zq1,2.,masse,w,zqm1)
+             do l=2,nlayer
+              zq( l,iq)=zq1(l)
+              zqm(l,iq)=zqm1(l)
+             enddo
+            enddo
+
+c           Surface condensation affects low winds
+            if (zmflux(1).lt.0) then 
+                zum(1)= zu(1) *  (w(1)/masse(1))
+                zvm(1)= zv(1) *  (w(1)/masse(1))
+                if (w(1).gt.masse(1)) then ! ensure numerical stability
+                  zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
+                  zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
+                end if
+            end if
+                    
+            ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
+            zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
+            zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
+            do iq=1,nq
+             zqm(nlayer+1,iq)= zq(nlayer,iq)
+            enddo
+
+#ifdef MESOSCALE
+!!!! AS: This part must be commented in the mesoscale model
+!!!! AS: ... to avoid instabilities.
+!!!! AS: you have to compile with -DMESOSCALE to do so
+#else  
+c           Tendencies on T, U, V, Q 
+c           """"""""""""""""""""""""
+            DO l=1,nlayer
+               IF(.not. co2clouds) THEN
+c             Tendencies on T
+                zdtsig(ig,l) = (1/masse(l)) *
+     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
+     &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
+     &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
+               ELSE
+                zdtsig(ig,l) = (1/masse(l)) *
+     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
+     &        - zmflux(l+1)*(ztm(l+1) - ztc(l)))               
+               ENDIF
+c D.BARDET: for diagnotics 
+                zmflux3D(ig,l)=zmflux(l)
+                ztm3D(ig,l)=ztm(l)
+                ztc3D(ig,l)=ztc(l)
+                
+                pdtc(ig,l) =  pdtc(ig,l) + zdtsig(ig,l)
+
+c             Tendencies on U
+                pduc(ig,l)   = (1/masse(l)) *
+     &        ( zmflux(l)*(zum(l) - zu(l))
+     &        - zmflux(l+1)*(zum(l+1) - zu(l)) )
+
+
+c             Tendencies on V
+                pdvc(ig,l)   = (1/masse(l)) *
+     &        ( zmflux(l)*(zvm(l) - zv(l))
+     &        - zmflux(l+1)*(zvm(l+1) - zv(l)) )
+
+            END DO
+
+#endif
+
+c           Tendencies on Q
+            do iq=1,nq
+!              if (noms(iq).eq.'co2') then 
+              if (iq.eq.ico2) then 
+c               SPECIAL Case when the tracer IS CO2 :
+                DO l=1,nlayer
+                  pdqc(ig,l,iq)= (1/masse(l)) *
+     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
+     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
+     &           + zcondicea(ig,l)*(zq(l,iq)-1.) )
+                END DO
+              else
+                DO l=1,nlayer
+                  pdqc(ig,l,iq)= (1/masse(l)) *
+     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
+     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
+     &           + zcondicea(ig,l)*zq(l,iq) )
+                END DO
+              end if
+            enddo
+
+       end if ! if (condsub)
+      END DO  ! loop on ig 
+
+c ***************************************************************
+c   CO2 snow / clouds scheme
+c ***************************************************************
+
+      call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
+     &        zcondicea,zcondices,zfallice,pemisurf)
+
+c ***************************************************************
+c Ecriture des diagnostiques
+c ***************************************************************
+
+c     DO l=1,nlayer
+c        DO ig=1,ngrid
+c         Taux de cond en kg.m-2.pa-1.s-1
+c          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
+c          Taux de cond en kg.m-3.s-1
+c          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
+c        END DO
+c     END DO
+c     call WRITEDIAGFI(ngrid,'tconda1',
+c    &'Taux de condensation CO2 atmospherique /Pa',
+c    & 'kg.m-2.Pa-1.s-1',3,tconda1)
+c     call WRITEDIAGFI(ngrid,'tconda2',
+c    &'Taux de condensation CO2 atmospherique /m',
+c    & 'kg.m-3.s-1',3,tconda2)
+
+! output falling co2 ice in 1st layer:
+!      call WRITEDIAGFI(ngrid,'fallice',
+!     &'Precipitation of co2 ice',
+!     & 'kg.m-2.s-1',2,zfallice(1,1))
+
+#ifndef MESOSCALE 
+! Extra special case for surface temperature tendency pdtsrfc:
+! we want to fix the south pole temperature to CO2 condensation temperature
+         if (caps.and.(obliquit.lt.27.)) then
+           ! check if last grid point is the south pole
+           if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
+           ! NB: Updated surface pressure, at grid point 'ngrid', is
+           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+!             write(*,*) "newcondens: South pole: latitude(ngrid)=",
+!     &                                           latitude(ngrid)
+             ztcondsol(ngrid)=
+     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
+     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
+           endif
+         endif
+#endif
+
+      END SUBROUTINE co2condens
+
+
+
+c *****************************************************************
+      SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
+c
+c    
+c     Operateur de moyenne inter-couche pour calcul de transport type
+c     Van-Leer " pseudo amont " dans la verticale
+c    q,w sont des arguments d'entree  pour le s-pg ....
+c    masse : masse de la couche Dp/g
+c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
+c    pente_max = 2 conseillee
+c
+c
+c   --------------------------------------------------------------------
+      IMPLICIT NONE
+
+c
+c
+c
+c   Arguments:
+c   ----------
+      integer nlayer
+      real masse(nlayer),pente_max
+      REAL q(nlayer),qm(nlayer+1)
+      REAL w(nlayer+1)
+c
+c      Local 
+c   ---------
+c
+      INTEGER l
+c
+      real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
+      real sigw, Mtot, MQtot
+      integer m 
+c     integer ismax,ismin 
+
+
+c    On oriente tout dans le sens de la pression 
+c     W > 0 WHEN DOWN !!!!!!!!!!!!!
+
+      do l=2,nlayer
+            dzqw(l)=q(l-1)-q(l)
+            adzqw(l)=abs(dzqw(l))
+      enddo
+
+      do l=2,nlayer-1
+            if(dzqw(l)*dzqw(l+1).gt.0.) then
+                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
+            else
+                dzq(l)=0.
+            endif
+            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
+            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
+      enddo
+
+         dzq(1)=0.
+         dzq(nlayer)=0.
+
+       do l = 1,nlayer-1
+
+c         Regular scheme (transfered mass < layer mass)
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
+             sigw=w(l+1)/masse(l+1)
+             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
+          else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
+             sigw=w(l+1)/masse(l)
+             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
+
+c         Extended scheme (transfered mass > layer mass)
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          else if(w(l+1).gt.0.) then
+             m=l+1
+             Mtot = masse(m)
+             MQtot = masse(m)*q(m)
+             do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
+                m=m+1
+                Mtot = Mtot + masse(m)
+                MQtot = MQtot + masse(m)*q(m)
+             end do
+             if (m.lt.nlayer) then
+                sigw=(w(l+1)-Mtot)/masse(m+1)
+                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*
+     &          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
+             else
+                w(l+1) = Mtot
+                qm(l+1) = Mqtot / Mtot
+                write(*,*) 'top layer is disapearing !'
+                stop
+             end if
+          else      ! if(w(l+1).lt.0) 
+             m = l-1 
+             Mtot = masse(m+1)
+             MQtot = masse(m+1)*q(m+1)
+             if (m.gt.0) then ! because some compilers will have problems
+                              ! evaluating masse(0)
+              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
+                m=m-1
+                Mtot = Mtot + masse(m+1)
+                MQtot = MQtot + masse(m+1)*q(m+1)
+                if (m.eq.0) exit
+              end do
+             endif
+             if (m.gt.0) then
+                sigw=(w(l+1)+Mtot)/masse(m)
+                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* 
+     &          (q(m)-0.5*(1.+sigw)*dzq(m))  )
+             else
+                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
+             end if   
+          end if
+       enddo
+
+c     boundary conditions (not used in newcondens !!)
+c         qm(nlayer+1)=0.
+c         if(w(1).gt.0.) then
+c            qm(1)=q(1)
+c         else 
+c           qm(1)=0.
+c         end if
+
+      END SUBROUTINE vl1d
+      
+      END MODULE co2condens_mod
Index: trunk/LMDZ.MARS/libf/phymars/newcondens.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/newcondens.F	(revision 2008)
+++ 	(revision )
@@ -1,979 +1,0 @@
-      SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep,
-     $                  pcapcal,pplay,pplev,ptsrf,pt,
-     $                  pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
-     $                  piceco2,psolaralb,pemisurf,
-     $                  pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
-     $                  fluxsurf_sw,zls,
-     $                  zdqssed_co2,pcondicea_co2microp,
-     $                  zdtcloudco2)
-                                                   
-       use tracer_mod, only: noms
-       use surfdat_h, only: emissiv, phisfi
-       use geometry_mod, only: latitude ! grid point latitudes (rad)
-       use planete_h
-       USE comcstfi_h
-#ifndef MESOSCALE
-       USE vertical_layers_mod, ONLY: bp
-#endif
-       IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Condensation/sublimation of CO2 ice on the ground and in the
-c   atmosphere 
-c   (Scheme described in Forget et al., Icarus, 1998)
-c
-c   author:   Francois Forget     1994-1996
-c   ------
-c             adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
-c
-c   input:
-c   ----- 
-c    ngrid                 nombre de points de verticales
-c                          (toutes les boucles de la physique sont au
-c                          moins vectorisees sur ngrid)
-c    nlayer                nombre de couches
-c    pplay(ngrid,nlayer)   Pressure levels 
-c    pplev(ngrid,nlayer+1) Pressure levels 
-c    pt(ngrid,nlayer)      temperature (en K)
-c    ptsrf(ngrid)          temperature de surface
-c 
-c                    \
-c    pdt(ngrid,nlayer)\  derivee temporelle physique avant condensation
-c                     /  ou sublimation pour  pt,ptsrf
-c    pdtsrf(ngrid)   /
-c
-c   output:
-c   -------
-c
-c    pdpsrf(ngrid)      \  derivee temporelle physique (contribution de
-c    pdtc(ngrid,nlayer) /  la condensation ou sublimation) pour Ps,pt,ptsrf
-c    pdtsrfc(ngrid)    /
-c
-c   Entree/sortie
-c   -------------
-c    
-c    piceco2(ngrid) :      quantite de glace co2 au sol (kg/m2)
-c    psolaralb(ngrid,2) :  albedo au sol
-c    pemisurf(ngrid)     :  emissivite du sol             
-
-c
-c=======================================================================
-c
-c    0.  Declarations :
-c    ------------------
-c
-      include "callkeys.h"
-
-c-----------------------------------------------------------------------
-c    Arguments :
-c    ---------
-      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
-      INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
-      INTEGER,INTENT(IN) :: nq  ! number of tracers
-
-      REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
-      REAL,INTENT(IN) :: pcapcal(ngrid)
-      REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
-      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
-      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
-      REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
-                                           ! previous physical processes (K/s)
-      REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
-                                           ! from previous physical processes
-      REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
-                                           ! from previous physical processes
-      REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
-                                       ! previous physical processes (K/s)
-      REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
-      REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
-      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
-      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
-                                              ! previous physical processes
-
-      REAL,INTENT(IN) :: zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
-      REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
-      REAL,INTENT(IN) :: zdtcloudco2(ngrid,nlayer) ! tendency on temperature due to latent heat                 
-
-      REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
-      REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
-      REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
-      
-      ! tendencies due to CO2 condensation/sublimation:
-      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
-      REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
-      REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
-      REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
-      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
-      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
-      
-      ! added to calculate flux dependent albedo:
-      REAL,intent(in) :: fluxsurf_sw(ngrid,2)
-      real,intent(in) :: zls ! solar longitude (rad)
-
-c
-c    Local variables :
-c    -----------------
-
-c   variables used for albedo parametrization       
-c --------------------------------------------      
-      INTEGER i,j
-c      REAL Fluxmean(jjp1)
-      INTEGER l,ig,iq,icap,nmix
-      LOGICAL transparency, fluxdependent
-c   flag transparency if you want to make the co2ice semi-transparent       
-      PARAMETER(transparency=.true.)
-c   flag fluxdependent if you want the co2ice albedo to be dependent on
-c   the incident solar flux      
-      PARAMETER(fluxdependent=.false.)
-      REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid)
-      REAL zt(ngrid,nlayer)
-      REAL zcpi
-      REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
-      REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
-      REAL zdiceco2(ngrid)
-      REAL zcondicea(ngrid,nlayer)
-      REAL zcondices(ngrid)
-      REAL zfallice(ngrid,nlayer+1) , zfallheat
-      REAL zmflux(nlayer+1)
-      REAL zu(nlayer),zv(nlayer)
-      REAL zq(nlayer,nq),zq1(nlayer)
-      REAL ztsrf(ngrid)
-      REAL ztc(nlayer), ztm(nlayer+1) 
-      REAL zum(nlayer+1) , zvm(nlayer+1)
-      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
-      REAL masse(nlayer),w(nlayer+1)
-      REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
-      LOGICAL condsub(ngrid)
-
-      real :: emisref(ngrid)
-
-c variable speciale diagnostique
-      real tconda1(ngrid,nlayer)
-      real tconda2(ngrid,nlayer)
-c     REAL zdiceco2a(ngrid) ! for diagnostic only
-      real zdtsig (ngrid,nlayer)
-      real zdt (ngrid,nlayer)
-      real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
-! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
-! then condensation temperature is computed using partial pressure of CO2
-      logical,parameter :: improved_ztcond=.true.
-! Bound co2 (tracer) values...
-      logical,parameter :: bound_qco2=.false.
-      real,parameter :: qco2max=1.1
-      real,parameter :: qco2mini=0.1
-      real :: zqco2
-
-c   local saved variables
-      integer,save :: ico2 ! index of CO2 tracer
-      real,save :: qco2min,qco2,mmean
-      real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
-      real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
-      real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
-      REAL,SAVE :: acond,bcond,ccond
-!      REAL,SAVE :: albediceF(ngrid)
-      real,save :: m_co2, m_noco2, A , B
-
-      LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
-
-      integer flag
-
-c D.BARDET: to debug 
-      real ztc3D(ngrid,nlayer)
-      REAL ztm3D(ngrid,nlayer)
-      REAL zmflux3D(ngrid,nlayer)
-c----------------------------------------------------------------------
-
-c   Initialisation
-c   --------------
-c
-      ! AS: firstcall OK absolute
-      IF (firstcall) THEN
-         
-         bcond=1./tcond1mb
-         ccond=cpp/(g*latcond)
-         acond=r/latcond
-
-         firstcall=.false.
-         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
-         write(*,*) 'Newcondens: bound_qco2=',bound_qco2
-         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
-         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,nq
-             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
-c          minimum CO2 mix. ratio below which mixing occur with layer above: 
-           qco2min =0.75  
-         end if
-      ENDIF ! of IF (firstcall)
-      zcpi=1./cpp
-
-c
-c======================================================================
-c    Calcul of CO2 condensation sublimation 
-c    ============================================================
-c  
-c    Used variable :
-c       piceco2(ngrid)   :  amount of co2 ice on the ground (kg/m2)
-c       zcondicea(ngrid,l):  condensation rate in layer  l  (kg/m2/s)
-c       zcondices(ngrid):  condensation rate on the ground (kg/m2/s)
-c       zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
-c                           
-c       pdtc(ngrid,nlayer)  : dT/dt due to cond/sub
-c
-c
-c     Tendencies set to 0 (except pdtc)
-c     -------------------------------------
-      DO l=1,nlayer
-         DO ig=1,ngrid
-           zcondicea(ig,l) = 0.
-           zfallice(ig,l) = 0.
-           pduc(ig,l)  = 0
-           pdvc(ig,l)  = 0
-         END DO
-      END DO
-         
-      DO iq=1,nq
-      DO l=1,nlayer
-         DO ig=1,ngrid
-           pdqc(ig,l,iq)  = 0
-         END DO
-      END DO
-      END DO
-
-      DO ig=1,ngrid
-         zfallice(ig,nlayer+1) = 0.
-         zcondices(ig) = 0.
-         pdtsrfc(ig) = 0.
-         pdpsrf(ig) = 0.
-         condsub(ig) = .false.
-         zdiceco2(ig) = 0.
-      ENDDO
-      zfallheat=0
-
-c     *************************
-c     ATMOSPHERIC CONDENSATION
-c     *************************
-
-c     Compute CO2 Volume mixing ratio
-c     -------------------------------
-      if (improved_ztcond.and.(ico2.ne.0)) then
-         DO l=1,nlayer
-            DO ig=1,ngrid
-              qco2=pq(ig,l,ico2)+pdq(ig,l,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,l) = qco2*mmean/m_co2 
-            ENDDO
-         ENDDO
-      else
-         DO l=1,nlayer
-            DO ig=1,ngrid
-              vmr_co2(ig,l)=0.95
-            ENDDO
-         ENDDO
-      end if
-
-      IF (.NOT. co2clouds) then
-c     forecast of atmospheric temperature zt and frost temperature ztcond
-c     --------------------------------------------------------------------
-
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
-!            ztcond(ig,l)=1./(bcond-acond*log(.0095*pplay(ig,l)))
-            if (pplay(ig,l).ge.1e-4) then
-              ztcond(ig,l)=
-     &         1./(bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l)))
-            else
-              ztcond(ig,l)=0.0 !mars Monica
-            endif
-         ENDDO
-      ENDDO
-
-      ztcond(:,nlayer+1)=ztcond(:,nlayer)
- 
-c      Condensation/sublimation in the atmosphere
-c      ------------------------------------------
-c      (calcul of zcondicea , zfallice and pdtc)
-c
-      DO l=nlayer , 1, -1
-         DO ig=1,ngrid
-           pdtc(ig,l)=0.
-           IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
-               condsub(ig)=.true.
-               IF (zfallice(ig,l+1).gt.0) then  
-                 zfallheat=zfallice(ig,l+1)*
-     &           (pphi(ig,l+1)-pphi(ig,l) +
-     &          cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
-               ELSE
-                   zfallheat=0.
-               ENDIF
-               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
-               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
-     &                        *ccond*pdtc(ig,l)- zfallheat
-c              Case when the ice from above sublimes entirely
-c              """""""""""""""""""""""""""""""""""""""""""""""
-               IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
-                  pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
-     &              (ccond*(pplev(ig,l)-pplev(ig,l+1)))
-                  zcondicea(ig,l)= -zfallice(ig,l+1)
-               END IF
-
-               zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
-            END IF
-         ENDDO
-      ENDDO
-      
-      ELSE
-         
-        DO ig=1,ngrid 
-         zfallice(ig,1) = zdqssed_co2(ig)
-        ENDDO
-        DO l=nlayer , 1, -1
-            DO ig=1,ngrid 
-         zcondicea(ig,l) = pcondicea_co2microp(ig,l)* 
-     &                        (pplev(ig,l) - pplev(ig,l+1))/g
-            ENDDO
-        ENDDO
-        DO l=nlayer, 1, -1
-         DO ig=1, ngrid
-            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
-            pdtc(ig,l)=0.
-         ENDDO
-      ENDDO 
-      ENDIF ! if not co2clouds
-
-      call WRITEdiagfi(ngrid,"pdtc_atm",
-     &         "temperature tendency due to CO2 condensation",
-     &         " ",3,pdtc) 
-     
-       call WRITEdiagfi(ngrid,"zcondicea",
-     &         "",
-     &         " ",3,zcondicea)         
-      
-       call WRITEdiagfi(ngrid,"zfallice",
-     &         "",
-     &         " ",2,zfallice(ngrid,1))
-
-c     *************************
-c       SURFACE  CONDENSATION
-c     *************************
-
-c     forecast of ground temperature ztsrf and frost temperature ztcondsol
-c     --------------------------------------------------------------------
-      DO ig=1,ngrid
-         ztcondsol(ig)=
-     &          1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1)))
-         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
-      ENDDO
-
-c
-c      Condensation/sublimation on the ground 
-c      --------------------------------------
-c      (compute zcondices and pdtsrfc)
-c
-      DO ig=1,ngrid
-         IF(latitude(ig).lt.0) THEN
-           ! Southern hemisphere
-            icap=2
-         ELSE
-           ! Northern hemisphere
-            icap=1
-         ENDIF
-       
-c
-c        Loop on where we have  condensation/ sublimation 
-         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.   ! ground cond 
-     $       (zfallice(ig,1).NE.0.) .OR.           ! falling snow
-     $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
-     $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
-            condsub(ig) = .true. 
-
-            IF (zfallice(ig,1).gt.0) then  
-                 zfallheat=zfallice(ig,1)*
-     &           (pphi(ig,1)- phisfi(ig) +
-     &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
-            ELSE
-                 zfallheat=0.
-            ENDIF
-
-c           condensation or partial sublimation of CO2 ice
-c           """""""""""""""""""""""""""""""""""""""""""""""
-            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 
-     &      /(latcond*ptimestep)         - zfallheat
-            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-       
-c           If the entire CO_2 ice layer sublimes
-c           """"""""""""""""""""""""""""""""""""""""""""""""""""
-c           (including what has just condensed in the atmosphere)
-
-            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
-     &          -zcondices(ig))THEN
-               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
-               pdtsrfc(ig)=(latcond/pcapcal(ig))*
-     &         (zcondices(ig)+zfallheat)
-            END IF
-
-c           Changing CO2 ice amount and pressure :
-c           """"""""""""""""""""""""""""""""""""
-
-            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
-            piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
-            pdpsrf(ig) = -zdiceco2(ig)*g
-
-            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
-               PRINT*,'STOP in condens'
-               PRINT*,'condensing more than total mass'
-               PRINT*,'Grid point ',ig
-               PRINT*,'Ps = ',pplev(ig,1)
-               PRINT*,'d Ps = ',pdpsrf(ig)
-               STOP
-            ENDIF
-         END IF ! if there is condensation/sublimmation
-      ENDDO ! of DO ig=1,ngrid
-
-c  ********************************************************************
-c  Surface albedo and emissivity of the surface below the snow (emisref)
-c  ********************************************************************
-c     Prepare the case where albedo varies with insolation:
-c     ----------------------------------------------------
-!      if (fluxdependent) then
-!
-c       Calcul du flux moyen (zonal mean)
-!        do j=1,jjp1
-!           Fluxmean(j)=0      
-!           do i=1,iim
-!               ig=1+(j-2)*iim +i
-!               if(j.eq.1) ig=1
-!               if(j.eq.jjp1) ig=ngrid
-!               Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
-!     $            +fluxsurf_sw(ig,2)
-!           enddo
-!           Fluxmean(j)=Fluxmean(j)/float(iim)
-!        enddo
-!
-c       const A and B used to calculate the albedo which depends on solar flux 
-c       albedice=constA+constB*Flux
-c       constT = time step to calculate the solar flux when flux decreases         
-!         constA=0.26
-c        constA=0.33
-c        constA=0.186
-!         constB=0.00187
-!         constT=10
-!      endif ! of if (fluxdependent)
-
-!     Check that amont of CO2 ice is not problematic
-      DO ig=1,ngrid
-           if(.not.piceco2(ig).ge.0.) THEN
-              if(piceco2(ig).le.-5.e-8) print*,
-     $         'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
-               piceco2(ig)=0.
-           endif
-      ENDDO
-      
-!     Set albedo and emissivity of the surface
-!     ----------------------------------------
-      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
-
-c     Calcul de l'albedo
-c     ------------------
-!         do ig =1,ngrid
-!           IF(ig.GT.ngrid/2+1) THEN
-!              icap=2
-!           ELSE
-!              icap=1
-!           ENDIF
-!           IF(firstcall2) THEN
-!           albediceF(ig)=albedice(icap)
-!           ENDIF
-c       if there is still co2ice        ccccccccccccccccccccccc
-!           if (piceco2(ig).gt.0) then
-!              emisref(ig) = emisice(icap)
-
-c     if flux dependent albedo is used 
-c     --------------------------------
-!             if (fluxdependent) then
-!                j=INT((ig-2)/iim)+2
-!                if(ig.eq.1) j=1
-!                if(ig.eq.ngrid) j=jjp1
-c                albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
-c     $          constA+constB*250)
-!                  albediceF_new(ig)=constA+constB*Fluxmean(j)
-!                if (albediceF(ig).gt.albediceF_new(ig)) then
-!                   albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
-!     $             constT)*(albediceF_new(ig)-albediceF(ig))
-!                else
-!                   albediceF(ig)=albediceF_new(ig)
-!                endif
-c               if part of the ice is transparent
-c               slopy = pente de la droite:  alpha = y*co2ice/1620
-c               pour une valeur superieur a une epaisseur de glace donnee
-c               ici, epaisseur limite = 10cm
-!                if (transparency) then
-!                slopy=1/(1620*0.10)
-!                alpha=MIN(slopy*piceco2(ig),1.)
-!                psolaralb(ig,1) = alpha*albediceF(ig)
-!     $                           +(1-alpha)*albedodat(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                else
-!                psolaralb(ig,1) = albediceF(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                endif
-!             else
-c           transparency set to true and fluxdependent set to false            
-!                if (transparency) then
-!                slopy=1/(1620*0.10)
-!                alpha=MIN(slopy*piceco2(ig),1.)
-!                psolaralb(ig,1) = alpha*albedice(icap)
-!     $                           +(1-alpha)*albedodat(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                else
-c           simplest case: transparency and flux dependent set to false
-!                psolaralb(ig,1) = albedice(icap)
-!                psolaralb(ig,2) = albedice(icap)
-!                endif
-!             endif    
-c         no more co2ice, albedo = ground albedo
-!           else
-!              psolaralb(ig,1) = albedodat(ig)
-!              psolaralb(ig,2) = albedodat(ig)
-!              emisref(ig) = emissiv
-!              pemisurf(ig) = emissiv
-!           endif
-!         end do ! end of the ig loop
-
-! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
-      DO ig=1,ngrid
-        if (piceco2(ig).eq.0) then
-          pemisurf(ig)=emissiv
-        endif
-      ENDDO
-
-!         firstcall2=.false.
-c ***************************************************************
-c  Correction to account for redistribution between sigma or hybrid 
-c  layers when changing surface pressure (and warming/cooling
-c  of the CO2 currently changing phase).
-c *************************************************************
-
-      DO ig=1,ngrid
-        if (condsub(ig)) then
-           do l=1,nlayer
-             ztc(l)  =zt(ig,l)   +pdtc(ig,l)  *ptimestep
-             zu(l)   =pu(ig,l)   +pdu( ig,l)  *ptimestep
-             zv(l)   =pv(ig,l)   +pdv( ig,l)  *ptimestep
-            do iq=1,nq
-             zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
-            enddo
-           end do
-
-c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
-c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-
-            zmflux(1) = -zcondices(ig)
-            DO l=1,nlayer
-             zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
-#ifndef MESOSCALE
-     &        + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
-c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
-          if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
-#else
-          if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
-#endif
-            END DO
-
-c Mass of each layer
-c ------------------ 
-            DO l=1,nlayer
-              masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
-            END DO
-
-
-c  Corresponding fluxes for T,U,V,Q
-c  """"""""""""""""""""""""""""""""
-
-c           averaging operator for TRANSPORT  
-c           """"""""""""""""""""""""""""""""
-c           Value transfert at the surface interface when condensation
-c           sublimation:
-            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
-            zum(1) = 0  
-            zvm(1) = 0  
-            do iq=1,nq
-              zqm(1,iq)=0. ! most tracer do not condense !
-            enddo
-c           Special case if one of the tracer is CO2 gas
-            if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
-
-c           Van Leer scheme:
-            DO l=1,nlayer+1
-                w(l)=-zmflux(l)*ptimestep
-            END DO
-            call vl1d(nlayer,ztc,2.,masse,w,ztm) 
-            call vl1d(nlayer,zu ,2.,masse,w,zum) 
-            call vl1d(nlayer,zv ,2.,masse,w,zvm) 
-            do iq=1,nq
-             do l=1,nlayer
-              zq1(l)=zq(l,iq)
-             enddo
-             zqm1(1)=zqm(1,iq)
-             call vl1d(nlayer,zq1,2.,masse,w,zqm1)
-             do l=2,nlayer
-              zq( l,iq)=zq1(l)
-              zqm(l,iq)=zqm1(l)
-             enddo
-            enddo
-
-c           Surface condensation affects low winds
-            if (zmflux(1).lt.0) then 
-                zum(1)= zu(1) *  (w(1)/masse(1))
-                zvm(1)= zv(1) *  (w(1)/masse(1))
-                if (w(1).gt.masse(1)) then ! ensure numerical stability
-                  zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
-                  zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
-                end if
-            end if
-                    
-            ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
-            zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
-            zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
-            do iq=1,nq
-             zqm(nlayer+1,iq)= zq(nlayer,iq)
-            enddo
-
-#ifdef MESOSCALE
-!!!! AS: This part must be commented in the mesoscale model
-!!!! AS: ... to avoid instabilities.
-!!!! AS: you have to compile with -DMESOSCALE to do so
-#else  
-c           Tendencies on T, U, V, Q 
-c           """"""""""""""""""""""""
-            DO l=1,nlayer
-               IF(.not. co2clouds) THEN
-c             Tendencies on T
-                zdtsig(ig,l) = (1/masse(l)) *
-     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
-     &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
-     &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
-               ELSE
-                zdtsig(ig,l) = (1/masse(l)) *
-     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
-     &        - zmflux(l+1)*(ztm(l+1) - ztc(l)))               
-               ENDIF
-c D.BARDET: for diagnotics 
-                zmflux3D(ig,l)=zmflux(l)
-                ztm3D(ig,l)=ztm(l)
-                ztc3D(ig,l)=ztc(l)
-                
-                pdtc(ig,l) =  pdtc(ig,l) + zdtsig(ig,l)
-
-c             Tendencies on U
-                pduc(ig,l)   = (1/masse(l)) *
-     &        ( zmflux(l)*(zum(l) - zu(l))
-     &        - zmflux(l+1)*(zum(l+1) - zu(l)) )
-
-
-c             Tendencies on V
-                pdvc(ig,l)   = (1/masse(l)) *
-     &        ( zmflux(l)*(zvm(l) - zv(l))
-     &        - zmflux(l+1)*(zvm(l+1) - zv(l)) )
-
-            END DO
-
-#endif
-
-c           Tendencies on Q
-            do iq=1,nq
-!              if (noms(iq).eq.'co2') then 
-              if (iq.eq.ico2) then 
-c               SPECIAL Case when the tracer IS CO2 :
-                DO l=1,nlayer
-                  pdqc(ig,l,iq)= (1/masse(l)) *
-     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
-     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
-     &           + zcondicea(ig,l)*(zq(l,iq)-1.) )
-                END DO
-              else
-                DO l=1,nlayer
-                  pdqc(ig,l,iq)= (1/masse(l)) *
-     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
-     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
-     &           + zcondicea(ig,l)*zq(l,iq) )
-                END DO
-              end if
-            enddo
-
-c           --------------------------------------------------------
-c           Roughly Simulate Molecular mixing when CO2 is too depleted by
-c           Surface condensation (mixing starts if qco2 < qco2min )  
-c           FF 06/2004
-c           WARNING : this is now done in convadj, better (FF 02/2005)
-c           --------------------------------------------------------
-            flag=0  ! now done in convadj : must be =0
-            if (flag.eq.1) then 
-            if(ico2.gt.0) then   ! relevant only if one tracer is CO2
-             if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
-     &       .lt.qco2min) then 
-                do iq=1,nq
-                  zq(1,iq)=pq(ig,1,iq)
-     &                     +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
-                  Smq(1,iq) = masse(1)*zq(1,iq)
-                end do
-                Sm(1)  = masse(1)
-                do l =2,nlayer
-                  do iq=1,nq
-                     zq(l,iq)=pq(ig,l,iq)
-     &                        +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
-                     smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
-                  end do 
-                  sm(l) = sm(l-1) + masse(l)
-                  if(zq(l,ico2).gt.qco2min) then 
-c                   mixmas: mass of atmosphere that must be mixed to reach qco2min
-                    mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
-     &                      /(zq(l,ico2)-qco2min)
-                    if((mixmas.le.sm(l)))then
-c                      OK if mixed mass less than mass of layers affected
-                       nmix=l   ! number of layer affected by mixing
-                       goto 99 
-                    end if
-                  end if
-                 end do
- 99              continue 
-                 do iq=1,nq
-                   qmix=zq(nmix,iq) 
-     &             +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
-                   do l=1,nmix-1
-                      pdqc(ig,l,iq)=
-     &                  (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
-                   end do 
-c                  layer only partly mixed :
-                   pdqc(ig,nmix,iq)=(
-     &             qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
-     &             -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
-
-                 end do 
-             end if
-            end if
-
-        endif   ! (flag.eq.1)
-       end if ! if (condsub)
-      END DO  ! loop on ig 
-
-c ***************************************************************
-c   CO2 snow / clouds scheme
-c ***************************************************************
-
-      call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
-     &        zcondicea,zcondices,zfallice,pemisurf)
-
-c ***************************************************************
-c Ecriture des diagnostiques
-c ***************************************************************
-
-c     DO l=1,nlayer
-c        DO ig=1,ngrid
-c         Taux de cond en kg.m-2.pa-1.s-1
-c          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
-c          Taux de cond en kg.m-3.s-1
-c          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
-c        END DO
-c     END DO
-c     call WRITEDIAGFI(ngrid,'tconda1',
-c    &'Taux de condensation CO2 atmospherique /Pa',
-c    & 'kg.m-2.Pa-1.s-1',3,tconda1)
-c     call WRITEDIAGFI(ngrid,'tconda2',
-c    &'Taux de condensation CO2 atmospherique /m',
-c    & 'kg.m-3.s-1',3,tconda2)
-
-! output falling co2 ice in 1st layer:
-!      call WRITEDIAGFI(ngrid,'fallice',
-!     &'Precipitation of co2 ice',
-!     & 'kg.m-2.s-1',2,zfallice(1,1))
-
-!! Specific stuff to bound co2 tracer values ....
-      if (bound_qco2.and.(ico2.ne.0)) then
-        do ig=1,ngrid
-          do l=1,nlayer
-            zqco2=pq(ig,l,ico2)
-     &            +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
-            if (zqco2.gt.qco2max) then
-            ! correct pdqc:
-              pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
-     &                               -pdq(ig,l,ico2)
-              write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
-     &                   " so that co2 conc. does not exceed",qco2max
-              write(*,*) "   ig:",ig,"  l:",l
-            endif ! of if (zqco2.gt.qco2max)
-            if (zqco2.lt.qco2mini) then
-            ! correct pdqc:
-              pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
-     &                               -pdq(ig,l,ico2)
-              write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
-     &                   " so that co2 conc. is not less than",qco2mini
-              write(*,*) "   ig:",ig,"  l:",l
-            endif ! of if (zqco2.lt.qco2mini)
-          end do 
-        enddo
-      endif ! of if (bound_qco2.and.(ico2.ne.0)) then
-
-#ifndef MESOSCALE 
-! Extra special case for surface temperature tendency pdtsrfc:
-! we want to fix the south pole temperature to CO2 condensation temperature
-         if (caps.and.(obliquit.lt.27.)) then
-           ! check if last grid point is the south pole
-           if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
-           ! NB: Updated surface pressure, at grid point 'ngrid', is
-           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
-!             write(*,*) "newcondens: South pole: latitude(ngrid)=",
-!     &                                           latitude(ngrid)
-             ztcondsol(ngrid)=
-     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
-     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
-             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
-           endif
-         endif
-#endif
-
-      end
-
-
-
-c *****************************************************************
-      SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
-c
-c    
-c     Operateur de moyenne inter-couche pour calcul de transport type
-c     Van-Leer " pseudo amont " dans la verticale
-c    q,w sont des arguments d'entree  pour le s-pg ....
-c    masse : masse de la couche Dp/g
-c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
-c    pente_max = 2 conseillee
-c
-c
-c   --------------------------------------------------------------------
-      IMPLICIT NONE
-
-c
-c
-c
-c   Arguments:
-c   ----------
-      integer nlayer
-      real masse(nlayer),pente_max
-      REAL q(nlayer),qm(nlayer+1)
-      REAL w(nlayer+1)
-c
-c      Local 
-c   ---------
-c
-      INTEGER l
-c
-      real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
-      real sigw, Mtot, MQtot
-      integer m 
-c     integer ismax,ismin 
-
-
-c    On oriente tout dans le sens de la pression 
-c     W > 0 WHEN DOWN !!!!!!!!!!!!!
-
-      do l=2,nlayer
-            dzqw(l)=q(l-1)-q(l)
-            adzqw(l)=abs(dzqw(l))
-      enddo
-
-      do l=2,nlayer-1
-            if(dzqw(l)*dzqw(l+1).gt.0.) then
-                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
-            else
-                dzq(l)=0.
-            endif
-            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
-            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
-      enddo
-
-         dzq(1)=0.
-         dzq(nlayer)=0.
-
-       do l = 1,nlayer-1
-
-c         Regular scheme (transfered mass < layer mass)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
-             sigw=w(l+1)/masse(l+1)
-             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
-          else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
-             sigw=w(l+1)/masse(l)
-             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
-
-c         Extended scheme (transfered mass > layer mass)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          else if(w(l+1).gt.0.) then
-             m=l+1
-             Mtot = masse(m)
-             MQtot = masse(m)*q(m)
-             do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
-                m=m+1
-                Mtot = Mtot + masse(m)
-                MQtot = MQtot + masse(m)*q(m)
-             end do
-             if (m.lt.nlayer) then
-                sigw=(w(l+1)-Mtot)/masse(m+1)
-                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*
-     &          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
-             else
-                w(l+1) = Mtot
-                qm(l+1) = Mqtot / Mtot
-                write(*,*) 'top layer is disapearing !'
-                stop
-             end if
-          else      ! if(w(l+1).lt.0) 
-             m = l-1 
-             Mtot = masse(m+1)
-             MQtot = masse(m+1)*q(m+1)
-             if (m.gt.0) then ! because some compilers will have problems
-                              ! evaluating masse(0)
-              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
-                m=m-1
-                Mtot = Mtot + masse(m+1)
-                MQtot = MQtot + masse(m+1)*q(m+1)
-                if (m.eq.0) exit
-              end do
-             endif
-             if (m.gt.0) then
-                sigw=(w(l+1)+Mtot)/masse(m)
-                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* 
-     &          (q(m)-0.5*(1.+sigw)*dzq(m))  )
-             else
-                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
-             end if   
-          end if
-       enddo
-
-c     boundary conditions (not used in newcondens !!)
-c         qm(nlayer+1)=0.
-c         if(w(1).gt.0.) then
-c            qm(1)=q(1)
-c         else 
-c           qm(1)=0.
-c         end if
-
-      end
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2008)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2009)
@@ -16,4 +16,5 @@
       use watercloud_mod, only: watercloud
       use watersat_mod, only: watersat
+      use co2condens_mod, only: co2condens
       use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
      &                        mem_Nccn_co2
@@ -1391,5 +1392,5 @@
                pdt(1:ngrid,1:nlayer) = 
      &              pdt(1:ngrid,1:nlayer) + 
-     &              zdtcloudco2(1:ngrid,1:nlayer)! --> in newcondens
+     &              zdtcloudco2(1:ngrid,1:nlayer)! --> in co2condens
 c            endif
             
@@ -1689,5 +1690,5 @@
 
       IF (callcond) THEN
-         CALL newcondens(ngrid,nlayer,nq,ptimestep,
+         CALL co2condens(ngrid,nlayer,nq,ptimestep,
      $              capcal,zplay,zplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
@@ -1765,5 +1766,5 @@
       IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
 !#ifndef MESOSCALE 
-!         if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
+!         if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
            ! NB: Updated surface pressure, at grid point 'ngrid', is
            !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
@@ -1777,5 +1778,5 @@
 c       everywhere except on the north permanent cap and in regions
 c       covered by dry ice. 
-c              ALWAYS PLACE these lines after newcondens !!!
+c              ALWAYS PLACE these lines after co2condens !!!
 c       -------------------------------------------------------------
          do ig=1,ngrid
