Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1920)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1921)
@@ -20,4 +20,5 @@
       USE newsedim_mod, ONLY: newsedim
       USE datafile_mod, ONLY: datadir
+
       IMPLICIT NONE
 
@@ -308,6 +309,7 @@
         firstcall=.false.
       ENDIF                     ! of IF (firstcall)
-   
-c-----Initialization
+c ===========================================================================   
+c     Initialization
+c ===========================================================================
       dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
       beta=0.85
@@ -335,7 +337,7 @@
       enddo
 
-c-------------------------------------------------------------------
+c ==========================================================================
 c   0.  Representation of sub-grid water ice clouds
-c-------------------------------------------------------------------
+c ==========================================================================
       IF (CLFvaryingCO2) THEN
 
@@ -345,5 +347,5 @@
          co2cloudfrac(:,:)=mincloud
          
-c-----Tendencies
+c Tendencies
          DO l=1,nlay
             DO ig=1,ngrid
@@ -439,6 +441,6 @@
          ENDDO ! of DO ig=1,ngrid
 !     Totalcloud frac of the column missing here
-c-----------------------
-c-----No sub-grid cloud representation (CLFvarying=false)
+c
+c No sub-grid cloud representation (CLFvarying=false)
       ELSE
          DO l=1,nlay
@@ -448,5 +450,5 @@
          END DO
       END IF                    ! end if (CLFvaryingco2)
-c------------------------------------------------------------------
+c =============================================================================
 c microtimestep timeloop for microphysics:
 c 0.Stepped entry for tendancies
@@ -454,7 +456,7 @@
 c 2.Call co2clouds microphysics
 c 3.Update tendancies 
-c------------------------------------------------------------------
+c =============================================================================
       DO microstep=1,imicroco2
-c------ Temperature tendency subpdt
+c Temperature tendency subpdt
         ! If imicro=1 subpdt is the same as pdt
         DO l=1,nlay
@@ -482,8 +484,10 @@
      &              sum_subpdq(ig,l,igcm_co2)
      &              + pdq(ig,l,igcm_co2)
-
+c D.BARDET :
+               if (co2useh2o) then
                sum_subpdq(ig,l,igcm_h2o_ice) = 
      &              sum_subpdq(ig,l,igcm_h2o_ice)
      &              + pdq(ig,l,igcm_h2o_ice)
+
                sum_subpdq(ig,l,igcm_ccn_mass) = 
      &              sum_subpdq(ig,l,igcm_ccn_mass)
@@ -492,7 +496,8 @@
      &              sum_subpdq(ig,l,igcm_ccn_number)
      &              + pdq(ig,l,igcm_ccn_number)
+                endif
           ENDDO
         ENDDO
-c- Effective tracers quantities in the cloud fraction
+c Effective tracers quantities in the cloud fraction
         IF (CLFvaryingCO2) THEN     
             pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
@@ -507,8 +512,8 @@
         END IF  
        
-c------------------------------------------------------
+c ========================================================================
 c 1.SEDIMENTATION : update tracers, compute parameters,
 c   call to sedimentation routine, update tendancies
-c------------------------------------------------------
+c ========================================================================
         IF (sedimentation) THEN
         
@@ -543,26 +548,30 @@
           ENDDO
         ENDDO
-!     Gravitational sedimentation       
+! Gravitational sedimentation       
         zqsed0(:,:,igcm_co2_ice)=zqsed(:,:,igcm_co2_ice)
         zqsed0(:,:,igcm_ccnco2_mass)=zqsed(:,:,igcm_ccnco2_mass)
         zqsed0(:,:,igcm_ccnco2_number)=zqsed(:,:,igcm_ccnco2_number) 
-       !We save actualized tracer values to compute sedimentation tendancies
+! We save actualized tracer values to compute sedimentation tendancies
         call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,ztsed,
      &     rsedcloudco2,rhocloudco2t,
      &     zqsed(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
-!     sedim at the surface of co2 ice : keep track of it for physiq_mod
+
+! sedim at the surface of co2 ice : keep track of it for physiq_mod
         do ig=1,ngrid 
           sum_subpdqs_sedco2(ig)=
      &         sum_subpdqs_sedco2(ig)+ wq(ig,1)/microtimestep
         end do
+
         call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,ztsed,
      &     rsedcloudco2,rhocloudco2t,
      &     zqsed(:,:,igcm_ccnco2_mass),wq,beta) 
+
         call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,ztsed,
      &     rsedcloudco2,rhocloudco2t,
      &     zqsed(:,:,igcm_ccnco2_number),wq,beta) 
+
         DO l = 1, nlay            !Compute tendencies
           DO ig=1,ngrid
@@ -578,5 +587,5 @@
           ENDDO
         ENDDO
-      !update subtimestep tendencies with sedimentation input
+! update subtimestep tendencies with sedimentation input
         DO l=1,nlay
          DO ig=1,ngrid
@@ -595,14 +604,14 @@
         END IF !(end if sedimentation)
         
-c------------------------------------------------------
+c ==============================================================================
 c      2.  Main call to the cloud schemes:
-c------------------------------------------------------
+c ==============================================================================
         CALL improvedCO2clouds(ngrid,nlay,microtimestep,
      &     pplay,pplev,pteff,sum_subpdt,
      &     pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
      &     nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
-c-----------------------------------------------------
+c ==============================================================================
 c      3.  Updating tendencies after cloud scheme:
-c-----------------------------------------------------
+c ==============================================================================
         DO l=1,nlay
           DO ig=1,ngrid
@@ -630,8 +639,10 @@
      &              sum_subpdq(ig,l,igcm_co2)
      &              + subpdqcloudco2(ig,l,igcm_co2)
-
+c D.BARDET :
+               if (co2useh2o) then
                sum_subpdq(ig,l,igcm_h2o_ice) =
      &              sum_subpdq(ig,l,igcm_h2o_ice)
      &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
+
                sum_subpdq(ig,l,igcm_ccn_mass) =
      &              sum_subpdq(ig,l,igcm_ccn_mass)
@@ -640,4 +651,5 @@
      &              sum_subpdq(ig,l,igcm_ccn_number)
      &              + subpdqcloudco2(ig,l,igcm_ccn_number)
+                endif
           ENDDO
         ENDDO
@@ -651,5 +663,5 @@
          pdqs_sedco2(ig)=sum_subpdqs_sedco2(ig)/real(imicroco2)
       enddo
-c------ Temperature tendency pdtcloud
+c Temperature tendency pdtcloud
       DO l=1,nlay
         DO ig=1,ngrid
@@ -658,5 +670,5 @@
         ENDDO
       ENDDO
-c------ Tracers tendencies pdqcloud
+c Tracers tendencies pdqcloud
       DO l=1,nlay
         DO ig=1,ngrid        
@@ -667,30 +679,31 @@
      &            sum_subpdq(ig,l,igcm_co2)/real(imicroco2)
      &            - pdq(ig,l,igcm_co2)
+c D.BARDET :
+             if (co2useh2o) then
              pdqcloudco2(ig,l,igcm_h2o_ice) = 
      &            sum_subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
      &            - pdq(ig,l,igcm_h2o_ice)
-        ENDDO
-      ENDDO
-      DO l=1,nlay
-        DO ig=1,ngrid
+
+             pdqcloudco2(ig,l,igcm_ccn_mass) = 
+     &            sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_mass)
+
+             pdqcloudco2(ig,l,igcm_ccn_number) = 
+     &            sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_number)
+             endif
+        
              pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
      &            sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
      &            - pdq(ig,l,igcm_ccnco2_mass)
+        
              pdqcloudco2(ig,l,igcm_ccnco2_number) = 
      &            sum_subpdq(ig,l,igcm_ccnco2_number)/real(imicroco2)
      &            - pdq(ig,l,igcm_ccnco2_number)
-             pdqcloudco2(ig,l,igcm_ccn_mass) = 
-     &            sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
-     &            - pdq(ig,l,igcm_ccn_mass)
-             pdqcloudco2(ig,l,igcm_ccn_number) = 
-     &            sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
-     &            - pdq(ig,l,igcm_ccn_number)
-        ENDDO
-      ENDDO
-      DO l=1,nlay
-        DO ig=1,ngrid
+
              pdqcloudco2(ig,l,igcm_dust_mass) = 
      &            sum_subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
      &            - pdq(ig,l,igcm_dust_mass)
+
              pdqcloudco2(ig,l,igcm_dust_number) = 
      &            sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2)
@@ -698,6 +711,6 @@
         ENDDO
       ENDDO
-c-------Due to stepped entry, other processes tendencies can add up to negative values
-c-------Therefore, enforce positive values and conserve mass
+c Due to stepped entry, other processes tendencies can add up to negative values
+c Therefore, enforce positive values and conserve mass
       DO l=1,nlay
         DO ig=1,ngrid
@@ -713,9 +726,12 @@
      &               - pqeff(ig,l,igcm_ccnco2_number)/ptimestep 
      &               - pdq(ig,l,igcm_ccnco2_number)+1.
+                
                 pdqcloudco2(ig,l,igcm_dust_number) =  
      &               -pdqcloudco2(ig,l,igcm_ccnco2_number)
+
                 pdqcloudco2(ig,l,igcm_ccnco2_mass) =
      &               - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
      &               - pdq(ig,l,igcm_ccnco2_mass)+1.e-20
+
                 pdqcloudco2(ig,l,igcm_dust_mass) = 
      &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
@@ -735,9 +751,12 @@
      &               - pqeff(ig,l,igcm_dust_number)/ptimestep 
      &               - pdq(ig,l,igcm_dust_number)+1.
+
                 pdqcloudco2(ig,l,igcm_ccnco2_number) =  
      &               -pdqcloudco2(ig,l,igcm_dust_number)
+
                 pdqcloudco2(ig,l,igcm_dust_mass) =
      &               - pqeff(ig,l,igcm_dust_mass)/ptimestep
      &               - pdq(ig,l,igcm_dust_mass) +1.e-20
+
                 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
      &               -pdqcloudco2(ig,l,igcm_dust_mass)
@@ -745,5 +764,5 @@
         ENDDO
       ENDDO
-        !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq      
+! pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq      
       DO l=1,nlay
         DO ig=1,ngrid
@@ -812,23 +831,25 @@
               endif
      
-      !update rice water
-          call updaterice_micro(
-     &    pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
-     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
-     &    pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
-     &    pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
-     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
-     &    pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
-     &    pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
-     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
-     &    pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
-     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-          
+c D.BARDET : update rice water only if co2 use h2o ice as CCN
+          if (co2useh2o) then
+             call updaterice_micro(
+     &       pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &       (pdq(ig,l,igcm_h2o_ice) +                     ! ice mass
+     &       pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &       pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
+     &       (pdq(ig,l,igcm_ccn_mass) +                    ! ccn mass
+     &       pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
+     &       pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
+     &       (pdq(ig,l,igcm_ccn_number) +                  ! ccn number
+     &       pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
+     &       tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+          endif
+
           call updaterdust(
      &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
-     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
+     &   (pdq(ig,l,igcm_dust_mass) +                     ! dust mass
      &    pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
      &    pqeff(ig,l,igcm_dust_number) +                 ! dust number
-     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
+     &   (pdq(ig,l,igcm_dust_number) +                   ! dust number
      &    pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
      &    rdust(ig,l))
@@ -881,23 +902,36 @@
         DO l=1,nlay
           DO ig=1,ngrid
+
+            pdqcloudco2(ig,l,igcm_ccnco2_mass)=
+     &        pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_ccnco2_number)=
+     &        pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_dust_mass)=
+     &        pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_dust_number)=
+     &        pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
+c D.BARDET
+            if (co2useh2o) then
+            pdqcloudco2(ig,l,igcm_h2o_ice)=
+     &        pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
+
             pdqcloudco2(ig,l,igcm_ccn_mass)=
      &        pdqcloudco2(ig,l,igcm_ccn_mass)*co2cloudfrac(ig,l)
-            pdqcloudco2(ig,l,igcm_ccnco2_mass)=
-     &        pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
+
             pdqcloudco2(ig,l,igcm_ccn_number)=
-     &        pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)
-            pdqcloudco2(ig,l,igcm_ccnco2_number)=
-     &        pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
-            pdqcloudco2(ig,l,igcm_dust_mass)=
-     &        pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
-            pdqcloudco2(ig,l,igcm_dust_number)=
-     &        pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
-            pdqcloudco2(ig,l,igcm_h2o_ice)=
-     &        pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
+     &        pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)            
+            endif
+
             pdqcloudco2(ig,l,igcm_co2_ice)=
      &        pdqcloudco2(ig,l,igcm_co2_ice)*co2cloudfrac(ig,l)
+
             pdqcloudco2(ig,l,igcm_co2)=
      &        pdqcloudco2(ig,l,igcm_co2)*co2cloudfrac(ig,l)
+
             pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*co2cloudfrac(ig,l)
+
             Qext1bins2(ig,l)=Qext1bins2(ig,l)*co2cloudfrac(ig,l)
           ENDDO
@@ -926,4 +960,9 @@
       call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
      &        ," ",2,tau1mic)
-         
+      call WRITEDIAGFI(ngrid,"memdNNccn","Nombre de CCN de glace d eau"
+     &        ,"kg/kg ",3,memdNNccn) 
+      call WRITEDIAGFI(ngrid,"memdMMccn","Masse de CCN de glace d eau"
+     &        ,"kg/kg ",3,memdMMccn) 
+      call WRITEDIAGFI(ngrid,"memdMMh2o","Masse de CCN de glace d eau"
+     &        ,"kg/kg ",3,memdMMh2o)         
       END
Index: trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1920)
+++ trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1921)
@@ -83,6 +83,6 @@
       LOGICAL,SAVE :: firstcall=.true.
       REAL*8   derf ! Error function
-      INTEGER ig,l,i,flag_pourri
-      
+      INTEGER ig,l,i
+
       real masse (ngrid,nlay) ! Layer mass (kg.m-2)
       REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
@@ -102,5 +102,5 @@
       DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
       DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
-      DOUBLE PRECISION Mo,No
+      DOUBLE PRECISION Mo,No,No_dust,Mo_dust
       DOUBLE PRECISION  Rn, Rm, dev2,dev3, n_derf, m_derf
       DOUBLE PRECISION memdMMccn(ngrid,nlay)
@@ -111,5 +111,4 @@
       DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
       DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
-      DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
       DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
 
@@ -138,23 +137,23 @@
       DOUBLE PRECISION ratioh2o_ccn
       DOUBLE PRECISION vo2co2
-c     Parameters of the size discretization
-c       used by the microphysical scheme
-      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum bounary radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
-      DOUBLE PRECISION vrat_cld ! Volume ratio
-      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
+
+c     Parameters of the size discretization used by the microphysical scheme
+      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9   ! Minimum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6   ! Maximum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10  ! Minimum bounary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4  ! Maximum boundary radius (m)
+      DOUBLE PRECISION vrat_cld                         ! Volume ratio
+      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1)         ! boundary values of each rad_cldco2 bin (m)
       SAVE rb_cldco2
-      DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
-      DOUBLE PRECISION vol_cld(nbinco2_cld)  ! particle volume for each bin (m3)
+      DOUBLE PRECISION dr_cld(nbinco2_cld)              ! width of each rad_cldco2 bin (m)
+      DOUBLE PRECISION vol_cld(nbinco2_cld)             ! particle volume for each bin (m3)
 
       DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
-      REAL sigma_iceco2 ! Variance of the co2 ice and CCN distributions
+      REAL sigma_iceco2   ! Variance of the co2 ice and CCN distributions
       SAVE sigma_iceco2 
-      REAL sigma_ice ! Variance of the h2o ice and CCN distributions
+      REAL sigma_ice      ! Variance of the h2o ice and CCN distributions
       SAVE sigma_ice
       DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
-      integer,save :: coeffh2o  !will be =0 is co2useh2o=.false.
+
 !     Variables for the meteoritic flux:
       integer,parameter :: nbin_meteor=100
@@ -217,9 +216,9 @@
 c        mteta  = 0.952
 c        mtetaco2 = 0.952
-        write(*,*) 'co2_param contact parameter:', mtetaco2
+c        write(*,*) 'co2_param contact parameter:', mtetaco2
+
 c       Volume of a co2 molecule (m3)
-        vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
-c       below, vo1 will be recomputed using real rho ice co2 (T-dependant)
-        vo1co2=vo1 ! AJOUT JA 
+        vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
+
 c       Variance of the ice and CCN distributions
         sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
@@ -231,5 +230,4 @@
         write(*,*) 'Volume of a co2 molecule:', vo1co2
         
-        coeffh2o=0
         if (co2useh2o) then
            write(*,*)
@@ -237,6 +235,6 @@
            write(*,*) "This means water ice particles can "
            write(*,*) "serve as CCN for CO2 microphysics"
-           coeffh2o=1
         endif
+
         if (meteo_flux) then
            write(*,*)
@@ -255,5 +253,4 @@
            endif
 !used Variables
-!           open(newunit=uMeteor,file=trim(datadir)//
            open(unit=uMeteor,file=trim(datadir)//
      &          '/Meteo_flux_Plane.dat'
@@ -298,5 +295,5 @@
 ! 1. Initialisation
 !=============================================================
-      flag_pourri=0
+
       meteor_ccn(:,:,:)=0.
       rice(:,:) = 1.e-8
@@ -373,14 +370,11 @@
            IF ( satu .ge. 1 ) THEN ! if there is condensation
 
-c     call updaterccn(zq(ig,l,igcm_dust_mass),
-c     &             zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
-
-              !We do rdust computation "by hand" because we don't want to
-              ! change the mininumum rdust value in updaterad... It would
-              ! mess up various part of the GCM !
+c We do rdust computation "by hand" because we don't want to
+c change the mininumum rdust value in updaterad... It would
+c mess up various part of the GCM !
               
               rdust(ig,l)= zq(ig,l,igcm_dust_mass)
-     &             *0.75/pi/rho_dust
-     &             / zq(ig,l,igcm_dust_number)
+     &             *0.75/(pi*rho_dust
+     &             * zq(ig,l,igcm_dust_number))
               rdust(ig,l)= rdust(ig,l)**(1./3.)
               if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20 
@@ -391,7 +385,15 @@
               rdust(ig,l)=min(5.e-4,rdust(ig,l))
       
-c       Expand the dust moments into a binned distribution
+c Expand the dust moments into a binned distribution
+              
+              n_aer(:)=0.
+              m_aer(:)=0.
+
               Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
               No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
+
+              No_dust=No
+              Mo_dust=Mo
+
               Rn = rdust(ig,l)
               Rn = -dlog(Rn) 
@@ -399,5 +401,5 @@
               n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
               m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
-      
+
               do i = 1, nbinco2_cld
                  n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
@@ -405,36 +407,49 @@
                  n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
                  m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
-                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
-     &                meteor_ccn(ig,l,i) !Ajout meteor_ccn particles
+                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
                  m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
-                 m_aer2(i)=4./3.*pi*rho_dust
+              enddo
+
+c Ajout meteor_ccn particles aux particules de poussière background
+              if (meteo_flux) then
+                 do i = 1, nbinco2_cld
+                    n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i) 
+                    m_aer(i) = m_aer(i) + 4./3.*pi*rho_dust
      &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
      &                *rad_cldco2(i)
-              enddo
-              if (meteo_flux) m_aer(:)=m_aer(:)+m_aer2(:)
+                 enddo
+              endif
               
-              !Same but with h2o particles as CCN
-              call updaterice_micro(zq(ig,l,igcm_h2o_ice),
-     &             zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
-     &             tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-              Mo = zq(ig,l,igcm_h2o_ice) +
-     &             zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30!*tauscaling(ig)
-   !  &             + 1.e-30 !Total mass of H20 crystals,CCN included
-              No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
-              Rn = rice(ig,l)
-              Rn = -dlog(Rn) 
-              Rm = Rn - 3. * sigma_ice*sigma_ice  
-              n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
-              m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
-              do i = 1, nbinco2_cld
-                 n_aer_h2oice(i) = -0.5 * No * n_derf 
-                 m_aer_h2oice(i) = -0.5 * Mo * m_derf 
-                 n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
-                 m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
-                 n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf 
-                 m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 
-                 rad_h2oice(i) = rad_cldco2(i)
-              enddo
-              !Call to nucleation routine
+c Same but with h2o particles as CCN only if co2useh2o=.true.
+              
+              n_aer_h2oice(:)=0.
+              m_aer_h2oice(:)=0.
+
+              if (co2useh2o) then
+                call updaterice_micro(zq(ig,l,igcm_h2o_ice),
+     &               zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
+     &                 tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+                Mo = zq(ig,l,igcm_h2o_ice) +
+     &               zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30 
+                     ! Total mass of H20 crystals,CCN included
+                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
+                Rn = rice(ig,l)
+                Rn = -dlog(Rn) 
+                Rm = Rn - 3. * sigma_ice*sigma_ice  
+                n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
+                m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
+                do i = 1, nbinco2_cld
+                  n_aer_h2oice(i) = -0.5 * No * n_derf 
+                  m_aer_h2oice(i) = -0.5 * Mo * m_derf 
+                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
+                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
+                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf 
+                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 
+                  rad_h2oice(i) = rad_cldco2(i)
+                enddo
+              endif
+             
+
+! Call to nucleation routine
               call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
      &             ,n_aer,rate,n_aer_h2oice
@@ -445,14 +460,18 @@
               dMh2o = 0.
               do i = 1, nbinco2_cld
-                 Proba=1.0-dexp(-1.*microtimestep*rate(i))
-                 Probah2o=coeffh2o*
-     &            (1.0-dexp(-1.*microtimestep*rateh2o(i))) !if co2useh2o=.false., this is =0
-                 dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
-                 dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
+                 Proba    =1.0-dexp(-1.*microtimestep*rate(i))
                  dN       = dN + n_aer(i) * Proba
                  dM       = dM + m_aer(i) * Proba             
               enddo
-        ! dM  masse activée (kg) et dN nb particules par  kg d'air
-        ! Now increment CCN tracers and update dust tracers
+              if (co2useh2o) then
+                 do i = 1, nbinco2_cld
+                    Probah2o = 1.0-dexp(-1.*microtimestep*rateh2o(i))
+                    dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
+                    dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
+                 enddo
+              endif
+
+! dM  masse activée (kg) et dN nb particules par  kg d'air
+! Now increment CCN tracers and update dust tracers
               dNN= dN/tauscaling(ig)
               dMM= dM/tauscaling(ig)
@@ -466,5 +485,6 @@
               zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
 
-c     Update CCN for CO2 nucleating on H2O CCN :
+
+c Update CCN for CO2 nucleating on H2O CCN :
               ! Warning: must keep memory of it 
               if (co2useh2o) then
@@ -489,8 +509,9 @@
                 memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
                 memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
-             endif !of if co2useh2o
-           ENDIF                ! of is satu >1
-!=============================================================
-! 4. Ice growth: scheme for radius evolution
+             endif ! of if co2useh2o
+           ENDIF   ! of is satu >1
+
+!=============================================================
+! 5. Ice growth: scheme for radius evolution
 !=============================================================
 
@@ -505,5 +526,4 @@
 
               Ic_rice=0.
-              flag_pourri=0
               lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
      &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
@@ -512,12 +532,11 @@
               !specific heat of atm cpp = 744.5 J.kg-1.K-1
 
-ccccccc  call scheme of microphys. mass growth for CO2
+c call scheme of microphys. mass growth for CO2
               call massflowrateCO2(pplay(ig,l),zt(ig,l),
      &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice) 
-c     Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance ! 
+c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance ! 
              
               if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then 
                  Ic_rice=0.
-                 flag_pourri=1
                  subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
                  dMice=0
@@ -531,5 +550,5 @@
               dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
               dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
-!facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
+! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
 ! latent heat release       >0 if growth i.e. if dMice >0
               subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
@@ -541,24 +560,24 @@
 
 !=============================================================
-! 5. Dust cores releasing if no more co2 ice :
+! 6. Dust cores releasing if no more co2 ice :
 !=============================================================
 
               if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
- !On sublime tout
+! On sublime tout
                  if (co2useh2o) then
-                 if (memdMMccn(ig,l) .gt. 0) then
+                   if (memdMMccn(ig,l) .gt. 0) then
                     zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
      &                   +memdMMccn(ig,l)
-                 endif
-                 if (memdMMh2o(ig,l) .gt. 0) then
+                   endif
+                   if (memdMMh2o(ig,l) .gt. 0) then
                     zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
      &                   +memdMMh2o(ig,l)
-                 endif
+                   endif
                  
-                 if (memdNNccn(ig,l) .gt. 0) then
+                   if (memdNNccn(ig,l) .gt. 0) then
                     zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
      &                   +memdNNccn(ig,l)
+                   endif
                  endif
-              endif
                     zq(ig,l,igcm_dust_mass) = 
      &                   zq(ig,l,igcm_dust_mass)
@@ -582,10 +601,10 @@
               endif !of if co2_ice <1e-25
 
-              ENDIF                ! of if NCCN > 1
+              ENDIF            ! of if NCCN > 1
           ENDDO                ! of ig loop
-        ENDDO                   ! of nlayer loop  
-
-!=============================================================
-! 6. END: get cloud tendencies 
+        ENDDO                  ! of nlayer loop  
+
+!=============================================================
+! 7. END: get cloud tendencies 
 !=============================================================
 
@@ -597,13 +616,17 @@
      &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
      &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
-     &       (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
-     &       zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
-     &       (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
-     &       zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
-     &       (zq(1:ngrid,1:nlay,igcm_ccn_number) -
-     &       zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
+        
+        if (co2useh2o) then
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
+     &         (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
+     &         zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
+     &         (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
+     &         zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
+     &         (zq(1:ngrid,1:nlay,igcm_ccn_number) -
+     &         zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
+        endif
+
         subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
      &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
@@ -619,4 +642,5 @@
      &       zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
 
+
         end
       
Index: trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F	(revision 1920)
+++ trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F	(revision 1921)
@@ -26,5 +26,7 @@
  ! nucrate_h2o en sortie aussi : 
 !nucleation sur dust et h2o separement ici
-#include "microphys.h"
+
+      include "microphys.h"
+      include "callkeys.h"
 
 c     Inputs
@@ -53,5 +55,5 @@
 
 
-      if (sat .gt. 1.) then    ! minimum condition to activate nucleation
+      IF (sat .gt. 1.) THEN    ! minimum condition to activate nucleation
 
         nco2   = pco2 / kbz / temp
@@ -63,5 +65,5 @@
 
 c       Loop over size bins
-        do 200 i=1,nbinco2_cld
+        do i=1,nbinco2_cld
 c            write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i),
 c     &          n_ccn(i)
@@ -70,23 +72,23 @@
 c           no dust, no need to compute nucleation!
             nucrate(i)=0.
-            goto 210
-          endif
+c            goto 210
+c          endif
+          else
+            if (rad_cldco2(i).gt.3000.*rstar) then
+              zefshapeco2 = fshapeco2simple
+            else
+             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
+            endif
 
-          if (rad_cldco2(i).gt.3000.*rstar) then
-            zefshapeco2 = fshapeco2simple
-          else
-            zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
-          endif
+            fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
+     &             zefshapeco2
+            deltaf = (2.*desorpco2-surfdifco2-fistar)/
+     &             (kbz*temp)
+            deltaf = min( max(deltaf, -100.d0), 100.d0)
 
-          fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
-     &             zefshapeco2
-          deltaf = (2.*desorpco2-surfdifco2-fistar)/
-     &             (kbz*temp)
-          deltaf = min( max(deltaf, -100.d0), 100.d0)
-
-          if (deltaf.eq.-100.) then
-            nucrate(i) = 0.
-          else
-            nucrate(i)= dble(sqrt ( fistar /
+            if (deltaf.eq.-100.) then
+                nucrate(i) = 0.
+            else
+                nucrate(i)= dble(sqrt ( fistar /
      &               (3.*pi*kbz*temp*(gstar*gstar)) )
      &                  * kbz * temp * rstar
@@ -98,31 +100,31 @@
 
             
-          endif
+            endif
+          endif ! if n_ccn(i) .lt. 1e-10
 
-210     continue
+          if (co2useh2o) then
 
-          if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
-c           no dust, no need to compute nucleation!
-            nucrate_h2oice(i)=0.
-            goto 200
-          endif
+            if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
+c           no H2O ice, no need to compute nucleation!
+               nucrate_h2oice(i)=0.
+            else   
+               if (rad_h2oice(i).gt.3000.*rstar) then
+                 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
+     &                (1.-mtetalocalh) / 4.
+               else  ! same m for dust/h2o ice
+                 zefshapeco2 = fshapeco2(mtetalocalh,
+     &                            (rad_h2oice(i)/rstar))
+               endif
 
-          if (rad_h2oice(i).gt.3000.*rstar) then
-            zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
-     &                (1.-mtetalocalh) / 4.
-          else
-            zefshapeco2 = fshapeco2(mtetalocalh,rad_h2oice(i)/rstar) ! same m for dust/h2o ice
-          endif
+               fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
+     &             zefshapeco2
+              deltaf = (2.*desorpco2-surfdifco2-fistar)/
+     &             (kbz*temp)
+              deltaf = min( max(deltaf, -100.d0), 100.d0)
 
-          fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
-     &             zefshapeco2
-          deltaf = (2.*desorpco2-surfdifco2-fistar)/
-     &             (kbz*temp)
-          deltaf = min( max(deltaf, -100.d0), 100.d0)
-
-          if (deltaf.eq.-100.) then
-            nucrate_h2oice(i) = 0.
-          else
-            nucrate_h2oice(i)= dble(sqrt ( fistar /
+              if (deltaf.eq.-100.) then
+                  nucrate_h2oice(i) = 0.
+              else
+                  nucrate_h2oice(i)= dble(sqrt ( fistar /
      &               (3.*pi*kbz*temp*(gstar*gstar)) )
      &                  * kbz * temp * rstar
@@ -132,11 +134,10 @@
      &                  / ( zefshapeco2 * nusco2 * m0co2 )
      &                  * dexp (deltaf))
+              endif
+            endif
           endif
-          
-         
+        enddo
 
-200     continue
-
-      else
+      ELSE ! parcelle d'air non saturée
 
         do i=1,nbinco2_cld
@@ -145,7 +146,6 @@
         enddo
 
-      endif
+      ENDIF ! if (sat .gt. 1.) 
 
-      return
       end
 
@@ -175,4 +175,3 @@
           fshapeco2 = 0.5*fshapeco2
 
-      return 
       end
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1920)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1921)
@@ -335,5 +335,5 @@
       REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
       REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
-      REAL icetotco2(ngrid)        ! Total mass of co2 ice (kg/m2) 
+      REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2) 
       REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
       REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
@@ -1284,49 +1284,49 @@
      &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
                 where (pq(:,:,igcm_ccn_mass) + 
-     &               ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
-                pdq(:,:,igcm_ccn_mass) = 
+     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
+                  pdq(:,:,igcm_ccn_mass) = 
      &               - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
-                pdq(:,:,igcm_ccn_number) = 
+                  pdq(:,:,igcm_ccn_number) = 
      &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
-             end where
-             where (pq(:,:,igcm_ccn_number) + 
-     &            ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
-             pdq(:,:,igcm_ccn_mass) = 
-     &            - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
-             pdq(:,:,igcm_ccn_number) = 
-     &            - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
-          end where
-       endif
+                end where
+                where (pq(:,:,igcm_ccn_number) + 
+     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
+                  pdq(:,:,igcm_ccn_mass) = 
+     &              - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+                  pdq(:,:,igcm_ccn_number) = 
+     &              - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+                end where
+             endif ! of if (co2useh2o)
 c     Negative values?
              where (pq(:,:,igcm_ccnco2_mass) + 
-     &            ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
-             pdq(:,:,igcm_ccnco2_mass) = 
+     &              ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
+                pdq(:,:,igcm_ccnco2_mass) = 
      &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
-             pdq(:,:,igcm_ccnco2_number) = 
+                pdq(:,:,igcm_ccnco2_number) = 
      &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
              end where
              where (pq(:,:,igcm_ccnco2_number) + 
-     &            ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
-             pdq(:,:,igcm_ccnco2_mass) = 
+     &              ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
+                pdq(:,:,igcm_ccnco2_mass) = 
      &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
-             pdq(:,:,igcm_ccnco2_number) = 
+                pdq(:,:,igcm_ccnco2_number) = 
      &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
-          end where
+             end where
        
 c     Negative values?
-          where (pq(:,:,igcm_dust_mass) + 
-     &         ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
-          pdq(:,:,igcm_dust_mass) = 
-     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
-          pdq(:,:,igcm_dust_number) = 
-     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
-       end where
-       where (pq(:,:,igcm_dust_number) + 
-     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
-       pdq(:,:,igcm_dust_mass) = 
-     &      - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
-       pdq(:,:,igcm_dust_number) = 
-     &      - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
-      end where
+             where (pq(:,:,igcm_dust_mass) + 
+     &              ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
+                pdq(:,:,igcm_dust_mass) = 
+     &            - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+                pdq(:,:,igcm_dust_number) = 
+     &            - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+             end where
+             where (pq(:,:,igcm_dust_number) + 
+     &              ptimestep*pdq(:,:,igcm_dust_number) < 0.)
+                pdq(:,:,igcm_dust_mass) = 
+     &            - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+                pdq(:,:,igcm_dust_number) = 
+     &            - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+             end where
       
       END IF                    ! of IF (co2clouds)
@@ -1367,6 +1367,6 @@
            zdqssed(1:ngrid,1:nq)=0
 
-cSedimentation for co2 clouds tracers are inside co2cloud microtimestep
-cZdqssed isn't
+c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
+c Zdqssed isn't
            call callsedim(ngrid,nlayer, ptimestep,
      &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
@@ -1374,5 +1374,5 @@
      &                pq, pdq, zdqsed, zdqssed,nq, 
      &                tau,tauscaling)
-   !Flux at the surface of co2 ice computed in co2cloud microtimestep
+c Flux at the surface of co2 ice computed in co2cloud microtimestep
            IF (co2clouds) THEN
               zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
@@ -1874,5 +1874,5 @@
                 enddo
               endif
-           endif
+           endif ! of (doubleq)
            if (co2clouds) then
               do ig=1,ngrid 
@@ -1882,17 +1882,19 @@
      &                zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
               enddo
+c D. BARDET compute integrated CO2 vapor and ice content
+              mtotco2(:)=0
+              icetotco2(:)=0
+              do ig=1,ngrid 
+               do l=1,nlayer
+                 mtotco2(ig) = mtotco2(ig) + 
+     &                      zq(ig,l,igcm_co2) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                 icetotco2(ig) = icetotco2(ig) + 
+     &                        zq(ig,l,igcm_co2_ice) * 
+     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+               enddo
+              enddo
               
-              
-              call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
-     &             'kg/kg',3,qccnco2)
-              call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
-     &             'part/kg',3,nccnco2)
-              call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
-     &             3,pq(:,:,igcm_co2_ice))
-              call WRITEDIAGFI(ngrid,'co2','co2','kg/kg',
-     &             3,pq(:,:,igcm_co2))
-                call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
-     &             3,pq(:,:,igcm_h2o_ice))
-           endif 
+           endif ! of if (co2clouds)
            
                     
@@ -2373,15 +2375,23 @@
           call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
      &                     "kg/kg",3,zq(1,1,igcm_co2))
-        
-         ! Compute co2 column
-         co2col(:)=0
-         do l=1,nlayer
-           do ig=1,ngrid
-             co2col(ig)=co2col(ig)+
-     &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
-           enddo
-         enddo
-         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
-     &                  co2col)
+          if (co2clouds) then
+            CALL WRITEDIAGFI(ngrid,'mtotco2',
+     &                       'total mass of CO2 vapor',
+     &                       'kg/m2',2,mtotco2)
+            CALL WRITEDIAGFI(ngrid,'zdtcloudco2',
+     &                       'temperature variation of CO2 latent heat',
+     &                       'K/s',3,zdtcloudco2)
+
+            CALL WRITEDIAGFI(ngrid,'icetotco2',
+     &                       'total mass of CO2 ice',
+     &                       'kg/m2',2,icetotco2)
+
+            call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
+     &                       'kg/kg',3,qccnco2)
+            call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
+     &                       'part/kg',3,nccnco2)
+            call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
+     &                       3,zq(:,:,igcm_co2_ice))
+          endif ! of if (co2clouds)
          endif ! of if (tracer.and.(igcm_co2.ne.0))
 
@@ -2435,4 +2445,9 @@
      &                       'Mean reff',
      &                       'm',2,rave)
+            call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
+     &             3,zq(:,:,igcm_h2o_ice))
+            call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
+     &             3,zq(:,:,igcm_h2o_vap))
+
 !A. Pottier
              CALL WRITEDIAGFI(ngrid,'rmoym',
