Index: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2447)
@@ -21,5 +21,5 @@
      &                      igcm_dust_submicron, rho_dust, rho_ice,
      &                      nqdust, igcm_stormdust_mass,
-     &                      igcm_topdust_mass
+     &                      igcm_topdust_mass, igcm_co2_ice
       use geometry_mod, only: latitude ! grid point latitudes (rad)
       use comgeomfi_h, only: sinlat ! sines of grid point latitudes
@@ -160,4 +160,6 @@
       CHARACTER(LEN=20) :: txt ! to temporarly store text
       CHARACTER(LEN=1) :: txt2 ! to temporarly store text
+! indexes of co2 ice :
+      INTEGER,SAVE :: i_co2ice=0 ! co2 ice
 ! indexes of dust scatterers:
       INTEGER,SAVE :: naerdust ! number of dust scatterers
@@ -194,4 +196,9 @@
           i_ice=igcm_h2o_ice
           write(*,*) "aeropacity: i_ice=",i_ice
+        ENDIF
+
+        IF (co2clouds.AND.activeco2ice) THEN
+          i_co2ice=igcm_co2_ice
+          write(*,*) "aeropacity: i_co2ice =",i_co2ice
         ENDIF
 
@@ -476,4 +483,46 @@
 
 c==================================================================
+        CASE("co2_ice") aerkind             ! CO2 ice crystals
+c==================================================================
+
+c       1. Initialization
+        aerosol(1:ngrid,1:nlayer,iaer) = 0.
+        taucloudvis(1:ngrid) = 0.
+        taucloudtes(1:ngrid) = 0.
+c       2. Opacity calculation
+        ! NO CLOUDS
+        IF (clearsky) THEN
+            aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
+        ! CLOUDSs
+        ELSE ! else (clearsky)
+          DO ig=1, ngrid
+            DO l=1,nlayer
+              aerosol(ig,l,iaer) = max(1E-20,
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_ice * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,i_co2ice) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+     &                              )
+              taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
+              taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
+     &          QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
+     &          ( 1.E0 - omegaREFir3d(ig,l,iaer) )
+            ENDDO
+          ENDDO
+          ! SUB-GRID SCALE CLOUDS
+          IF (CLFvarying) THEN
+             DO ig=1, ngrid
+                DO l=1,nlayer-1
+                   CLFtot  = max(totcloudfrac(ig),0.01) 
+                   aerosol(ig,l,iaer)=
+     &                    aerosol(ig,l,iaer)/CLFtot
+                   aerosol(ig,l,iaer) = 
+     &                    max(aerosol(ig,l,iaer),1.e-9)
+                ENDDO
+             ENDDO
+          ENDIF ! end (CLFvarying)             
+        ENDIF ! end (clearsky)
+
+c==================================================================
         CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for 
 c       stormdust  (transport of mass and number mixing ratio) 
Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2447)
@@ -14,6 +14,6 @@
      &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
-     &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
-     &   ,satindexco2,rdstorm,slpwind,calllott_nonoro                   &
+     &   ,co2useh2o,meteo_flux,activeco2ice,CLFvaryingCO2,spantCO2      &
+     &   ,CLFvarying,satindexco2,rdstorm,slpwind,calllott_nonoro        &
      &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
      &   ,hdo,hdofrac
@@ -72,4 +72,5 @@
       logical activice,tifeedback,supersat,caps
       logical co2clouds,co2useh2o,meteo_flux,CLFvaryingCO2,satindexco2
+      logical activeco2ice
       integer spantCO2
       logical CLFvarying
Index: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2447)
@@ -10,7 +10,7 @@
      $     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
      &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
-     $     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
-     $     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
-     $     clearsky,totcloudfrac)
+     $     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
+     $     rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
+     $     alpha_hmons,nohmons,clearsky,totcloudfrac)
 
       use aeropacity_mod, only: aeropacity
@@ -19,5 +19,5 @@
       use dimradmars_mod, only: naerkind, name_iaer,
      &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice,
+     &            iaer_dust_submicron, iaer_h2o_ice, iaer_co2_ice,
      &            iaer_stormdust_doubleq,iaer_topdust_doubleq
       use yomlw_h, only: gcp, nlaylte
@@ -203,4 +203,6 @@
       REAL,INTENT(OUT) :: rice(ngrid,nlayer)   ! Ice geometric mean radius (m)
       REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
+      REAL,INTENT(OUT) :: riceco2(ngrid,nlayer) ! CO2 ice mean radius(m)
+      REAL,INTENT(OUT) :: nuiceco2(ngrid,nlayer) ! Effective variance
       REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
 
@@ -317,4 +319,5 @@
          iaer_dust_submicron=0
          iaer_h2o_ice=0
+         iaer_co2_ice=0
          iaer_stormdust_doubleq=0
          iaer_topdust_doubleq=0
@@ -349,4 +352,12 @@
              if (name_iaer(iaer).eq."h2o_ice") then
                iaer_h2o_ice = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         endif
+         if (co2clouds.AND.activeco2ice) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."co2_ice") then
+               iaer_co2_ice = iaer
                aer_count = aer_count + 1
              endif
@@ -412,11 +423,9 @@
 c     Computing aerosol optical properties and opacity
 c     ------------------------------------------------
-
 c     Updating aerosol size distributions:
       CALL updatereffrad(ngrid,nlayer,
      &                rdust,rstormdust,rtopdust,rice,nuice,
-     &                reffrad,nueffrad,
+     &                reffrad,nueffrad, riceco2, nuiceco2,
      &                pq,tauscaling,tau,pplay)
-
 c     Computing 3D scattering parameters:
       CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
@@ -425,5 +434,4 @@
      &                      QREFvis3d,QREFir3d,
      &                      omegaREFvis3d,omegaREFir3d)
-
 c     Computing aerosol optical depth in each layer:
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
Index: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2447)
@@ -673,5 +673,4 @@
        ENDIF ! of IF(microphys)
       endif ! of if (water)
-
       END SUBROUTINE callsedim
       
Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F90	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F90	(revision 2447)
@@ -98,5 +98,5 @@
 !======================================================================================================================!
   subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
-                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rsedcloud, rhocloud, pzlev, pdqs_sedco2, pdu, &
+                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev, pdqs_sedco2, pdu, &
                       pu, pcondicea, co2ice)
 
@@ -115,4 +115,8 @@
 
   use improvedCO2clouds_mod, only: improvedCO2clouds, density_co2_ice
+
+#ifndef MESOSCALE
+  use vertical_layers_mod, only: ap, bp
+#endif
 
   implicit none
@@ -151,7 +155,6 @@
 !------------------
   real, intent(out) :: &
-     rdust(ngrid,nlay),         & ! Dust geometric mean radius (m)
      rice(ngrid,nlay),          & ! Water Ice mass mean radius (m)
-     rsedcloud(ngrid,nlay),     & ! Water Cloud sedimentation radius
+!     rsedcloud(ngrid,nlay),     & ! Water Cloud sedimentation radius
      rhocloud(ngrid,nlay),      & ! Water Cloud density (kg.m-3)
      pdqs_sedco2(ngrid),        & ! CO2 flux at the surface
@@ -159,5 +162,10 @@
      pcondicea(ngrid,nlay),     & ! Rate of condensation/sublimation of co2 ice in layers
      pdtcloudco2(ngrid,nlay),   & ! Tendency on temperature due to latent heat
-     riceco2(ngrid,nlay)          ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
+     riceco2(ngrid,nlay),       & ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
+     rsedcloudco2(ngrid,nlay)     ! Cloud sedimentation radius
+
+  real, intent(inout) :: &
+     rdust(ngrid,nlay) ! Dust geometric mean radius (m)
+
 !----------------------------------------------------------------------------------------------------------------------!
 ! Local:
@@ -225,5 +233,4 @@
      sum_subpdt(ngrid,nlay),        &! On temperature, cf. pdtcloud
 !  ---Sedimentation
-     rsedcloudco2(ngrid,nlay),      &! Cloud sedimentation radius
      ztsed(ngrid,nlay),             &! Temperature with real-time value in microtimeloop
      zqsed(ngrid,nlay,nq),          &! Tracers with real-time value in µloop
@@ -249,6 +256,10 @@
      wq(ngrid,nlay+1),              &! Displaced tracer mass (kg.m-2) during microtimestep
      satuco2(ngrid,nlay),           &! CO2 satu ratio for output diagfi
-     zqsatco2(ngrid,nlay)            ! Saturation co2
-
+     zqsatco2(ngrid,nlay),          &! Saturation co2
+     availco2,&
+     masslayer, &
+     tmp, a,b, &
+     new_pdq(ngrid,nlay)
+     
   double precision :: &
 ! ---Extinction coefficients at 1 micron of CO2 particles
@@ -416,5 +427,5 @@
   do l = 1, nlay
     do ig = 1, ngrid
-      masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
+      masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1)) ) / g
       epaisseur(ig,l) = pzlev(ig,l+1) - pzlev(ig,l)
     end do
@@ -545,47 +556,110 @@
 ! 4. Microphysics of CO2 cloud formation
 !----------------------------------------------------------------------------------------------------------------------!
+  pqeff(:,:,:) = pq(:,:,:)
+  pteff(:,:) = pt(:,:)
+!----------------------------------------------------------------------------------------------------------------------!
+! 4.2. Effective tracers quantities in the cloud
+!----------------------------------------------------------------------------------------------------------------------!
+  if (CLFvaryingCO2) then
+    pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:)
+
+    pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:)
+
+    pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:)
+  end if
+!----------------------------------------------------------------------------------------------------------------------!
   do microstep = 1, imicroco2
 !----------------------------------------------------------------------------------------------------------------------!
 ! 4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
+!----------------------------------------------------------------------------------------------------------------------!
+   do l = 1, nlay
+     do ig = 1, ngrid
+       ! on temperature
+       sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l)
+        
+       ! on tracers
+       sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)
+
+       sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)
+
+       sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass)
+
+       sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number)
+
+       sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice)
+
+       sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)
+
+       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) + pdq(ig,l,igcm_ccn_mass)
+
+         sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)
+       end if
+     end do ! ngrid
+   end do ! nlay
+!----------------------------------------------------------------------------------------------------------------------!
+! 4.4. Main call to the cloud scheme
+!----------------------------------------------------------------------------------------------------------------------!
+   call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
+                          subpdqcloudco2, subpdtcloudco2, nq, tauscaling, mem_Mccn_co2, mem_Mh2o_co2, mem_Nccn_co2, &
+                          rb_cldco2, sigma_iceco2, dev2)
+   do l = 1, nlay
+     do ig = 1, ngrid
+       if(pq(ig,l,igcm_co2_ice) + microtimestep*(sum_subpdq(ig,l,igcm_co2_ice)+subpdqcloudco2(ig,l,igcm_co2_ice)) &
+         .le. 1.e-12) then
+         subpdqcloudco2(ig,l,igcm_co2_ice) = -pq(ig,l,igcm_co2_ice)/microtimestep - sum_subpdq(ig,l,igcm_co2_ice)
+         subpdqcloudco2(ig,l,igcm_co2) = -subpdqcloudco2(ig,l,igcm_co2_ice)
+       end if
+
+       if(pq(ig,l,igcm_co2) + microtimestep*(sum_subpdq(ig,l,igcm_co2)+subpdqcloudco2(ig,l,igcm_co2)) .le. 1.e-12) then
+         subpdqcloudco2(ig,l,igcm_co2) =  - pq(ig,l,igcm_co2)/microtimestep - sum_subpdq(ig,l,igcm_co2)
+         subpdqcloudco2(ig,l,igcm_co2_ice) = -subpdqcloudco2(ig,l,igcm_co2)
+       end if
+  ! ccnco2_number and ccnco2_mass
+       if (((pq(ig,l,igcm_ccnco2_number)+(sum_subpdq(ig,l,igcm_ccnco2_number)+subpdqcloudco2(ig,l,igcm_ccnco2_number)) &
+          *microtimestep).le.1.) .or. &
+         (pq(ig,l,igcm_ccnco2_mass)+(sum_subpdq(ig,l,igcm_ccnco2_mass)+subpdqcloudco2(ig,l,igcm_ccnco2_mass)) &
+          *microtimestep.le.1e-20)) then
+         subpdqcloudco2(ig,l,igcm_ccnco2_number) = - pq(ig,l,igcm_ccnco2_number)/microtimestep + 1. &
+                                                  - sum_subpdq(ig,l,igcm_ccnco2_number)
+         subpdqcloudco2(ig,l,igcm_dust_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_number)
+
+         subpdqcloudco2(ig,l,igcm_ccnco2_mass) = - pq(ig,l,igcm_ccnco2_mass)/microtimestep + 1e-20 &
+                                                - sum_subpdq(ig,l,igcm_ccnco2_mass)
+         subpdqcloudco2(ig,l,igcm_dust_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_mass)
+       end if
+     end do
+   end do
+!----------------------------------------------------------------------------------------------------------------------!
+! 4.5. Updating tendencies after cloud scheme
 !----------------------------------------------------------------------------------------------------------------------!
     do l = 1, nlay
       do ig = 1, ngrid
-        ! on temperature
-        sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l)
-
-        ! on tracers
-        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)
-
-        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)
-
-        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass)
-
-        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number)
-
-        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice)
-
-        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)
+        sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
+
+        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
+
+        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
+
+        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
+
+        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
+                                              subpdqcloudco2(ig,l,igcm_ccnco2_number)
+
+        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
+
+        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
 
         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) + pdq(ig,l,igcm_ccn_mass)
-
-          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)
+          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) + subpdqcloudco2(ig,l,igcm_ccn_mass)
+
+          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
         end if
       end do ! ngrid
     end do ! nlay
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.2. Effective tracers quantities in the cloud
-!----------------------------------------------------------------------------------------------------------------------!
-    pqeff(:,:,:) = pq(:,:,:)
-
-    if (CLFvaryingCO2) then
-      pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:)
-
-      pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:)
-
-      pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:)
-    end if
 !----------------------------------------------------------------------------------------------------------------------!
 ! 4.3. Gravitational sedimentation
@@ -612,13 +686,10 @@
 
           ! Get density cloud and co2 ice particle radius
-          call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
-
-          ! conflict with updaterice_microco2 C.M.
-          if (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1) then
-            riceco2(ig,l) = 1.e-9
+          if (Niceco2.ne.0d0) then
+            call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
+          else
+            riceco2(ig,l) = 0.
+            rhocloudco2t(ig,l) = 0.
           end if
-
-          ! done in updaterice_microco2 C.M.
-          rhocloudco2t(ig,l) = min( max(rhocloudco2t(ig,l),rho_ice_co2T(ig,l)), rho_dust )
         end do ! ngrid
       end do ! nlay
@@ -632,5 +703,11 @@
 ! 4.3.c. Sedimentation of co2 ice
 !----------------------------------------------------------------------------------------------------------------------!
-      rsedcloudco2(1:ngrid, 1:nlay) = 5e-6
+      do ig = 1, ngrid
+        do l = 1, nlay
+          rsedcloudco2(ig,l) = max( riceco2(ig,l)*(1.+sigma_iceco2)*(1.+sigma_iceco2)*(1.+sigma_iceco2), &
+                                    rdust(ig,l) )
+        end do
+      end do
+
       wq(:,:) = 0.
       call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
@@ -638,5 +715,5 @@
 
       do ig = 1, ngrid
-        sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep
+        sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep !wq est en kg.m-2
       end do
 !----------------------------------------------------------------------------------------------------------------------!
@@ -653,5 +730,4 @@
       call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
                    rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_number), wq, beta)
-
       !TODO: ajouter le calcule de la tendance a la surface comme co2ice
 !----------------------------------------------------------------------------------------------------------------------!
@@ -669,5 +745,4 @@
         end do
       end do
-
       ! update subtimestep tendencies with sedimentation input
       do l = 1, nlay
@@ -681,44 +756,9 @@
       end do
     end if !(end if sedimentation)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.4. Main call to the cloud scheme
-!----------------------------------------------------------------------------------------------------------------------!
-    call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
-                           subpdqcloudco2, subpdtcloudco2, nq, tauscaling, mem_Mccn_co2, mem_Mh2o_co2, mem_Nccn_co2, &
-                           rb_cldco2, sigma_iceco2, dev2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.5. Updating tendencies after cloud scheme
-!----------------------------------------------------------------------------------------------------------------------!
-    do l = 1, nlay
-      do ig = 1, ngrid
-        sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
-
-        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
-
-        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
-
-        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
-
-        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
-                                              subpdqcloudco2(ig,l,igcm_ccnco2_number)
-
-        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
-
-        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
-
-        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) + subpdqcloudco2(ig,l,igcm_ccn_mass)
-
-          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
-        end if
-      end do ! ngrid
-    end do ! nlay
   end do ! of do microstep = 1, imicroco2
 !----------------------------------------------------------------------------------------------------------------------!
 ! 5. Compute final tendencies after time loop
 !----------------------------------------------------------------------------------------------------------------------!
-! condensation/sublimation rate in the atmosphere
+! condensation/sublimation rate in the atmosphere (kg.kg-1.s-1)
   do l = nlay, 1, -1
     do ig = 1, ngrid
@@ -732,5 +772,5 @@
   end do
 
-  ! temperature tendency
+  ! temperature tendency (T.s-1)
   do l = 1, nlay
     do ig = 1, ngrid
@@ -753,5 +793,5 @@
 
       pdqcloudco2(ig,l,igcm_dust_mass) = ( sum_subpdq(ig,l,igcm_dust_mass) / real(imicroco2) ) - &
-                                         pdq(ig,l,igcm_dust_mass)
+                                          pdq(ig,l,igcm_dust_mass)
 
       pdqcloudco2(ig,l,igcm_dust_number) = ( sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2) ) - &
@@ -759,9 +799,9 @@
 
       if (co2useh2o) then
-        pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) -&
-                                         pdq(ig,l,igcm_h2o_ice)
+        pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) - &
+                                          pdq(ig,l,igcm_h2o_ice)
 
         pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - &
-                                          pdq(ig,l,igcm_ccn_mass)
+                                           pdq(ig,l,igcm_ccn_mass)
 
         pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
@@ -778,4 +818,5 @@
     do ig = 1, ngrid
       Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep
+      Niceco2 = max (Niceco2, threshold)
 
       Nccnco2 = max( (pqeff(ig,l,igcm_ccnco2_number) + (pdq(ig,l,igcm_ccnco2_number) + &
@@ -848,12 +889,4 @@
   end do
 !----------------------------------------------------------------------------------------------------------------------!
-! 8. Compute water cloud sedimentation radius
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig= 1, ngrid
-      rsedcloud(ig,l) = max( rice(ig,l)*(1.+nuice_sed)* (1.+nuice_sed)*(1.+nuice_sed), rdust(ig,l) )
-    end do
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
 ! 9. CO2 saturated quantities
 !----------------------------------------------------------------------------------------------------------------------!
@@ -916,5 +949,7 @@
   call WRITEDIAGFI(ngrid, "satuco2", "vap in satu", " ", 3, satuco2)
 
-  call WRITEDIAGFI(ngrid, "co2cloudfrac", "co2 cloud fraction", " ", 3, co2cloudfrac)
+  if (CLFvaryingCO2) then
+    call WRITEDIAGFI(ngrid, "co2cloudfrac", "co2 cloud fraction", " ", 3, co2cloudfrac)
+  end if
 
   call writediagfi(ngrid, "riceco2", "ice radius", "m", 3, riceco2)
Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod4micro.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2condens_mod4micro.F90	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/co2condens_mod4micro.F90	(revision 2447)
@@ -194,4 +194,16 @@
   logical :: &
      condsub(ngrid) ! True if there is condensation/sublimation (used for co2snow)
+
+
+  ! check with co2 cloud parameterisation
+  real :: &
+     zt_2(ngrid,nlayer), &
+     ztcond_2(ngrid,nlayer+1), &
+     zfallice_2(ngrid,nlayer+1), &
+     pdtc_2(ngrid,nlayer), &
+     zfallheat_2, &
+     zcondicea_2(ngrid,nlayer), &
+     diff_zcondicea(ngrid,nlayer), &
+     diff_zfallice(ngrid)
 !======================================================================================================================!
 ! BEGIN ===============================================================================================================!
@@ -281,5 +293,5 @@
 ! Only sedimentation falls on the ground !
   do ig = 1, ngrid
-      zfallice(ig) =  zdqssed_co2(ig)
+      zfallice(ig) = zdqssed_co2(ig)
       piceco2(ig) = piceco2(ig) + zfallice(ig)*ptimestep
   end do
@@ -288,4 +300,49 @@
 
   call writediagfi(ngrid,"zfallice", "Sedimentation rate", " ", 2, zfallice)
+
+! Compute without microphysics
+ do l =1, nlayer
+   do ig = 1, ngrid
+     zt_2(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
+   end do
+ end do
+
+ do l = 1, nlayer
+   do ig = 1, ngrid
+     if (pplay(ig,l).ge.1e-4) then
+       ztcond_2(ig,l) = 1. / (bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l)))
+     else
+       ztcond_2(ig,l) = 0.0 !mars Monica
+     endif
+   end do
+ end do
+
+ ztcond_2(:,nlayer+1)=ztcond_2(:,nlayer)
+ zfallice_2(:,:) = 0
+
+ do l = nlayer , 1, -1
+   do ig = 1, ngrid
+     pdtc_2(ig,l) = 0.
+     if ((zt_2(ig,l).lt.ztcond_2(ig,l)).or.(zfallice_2(ig,l+1).gt.0)) then
+       if (zfallice_2(ig,l+1).gt.0) then  
+         zfallheat_2 = zfallice_2(ig,l+1)*(pphi(ig,l+1)-pphi(ig,l) + cpice*(ztcond_2(ig,l+1)-ztcond_2(ig,l)))/latcond
+       else
+         zfallheat_2 = 0.
+       end if
+       pdtc_2(ig,l) = (ztcond_2(ig,l) - zt_2(ig,l))/ptimestep
+       zcondicea_2(ig,l) = (pplev(ig,l)-pplev(ig,l+1))*ccond*pdtc(ig,l)- zfallheat_2
+       ! Case when the ice from above sublimes entirely
+       ! """""""""""""""""""""""""""""""""""""""""""""""
+       if (zfallice_2(ig,l+1).lt.- zcondicea_2(ig,l)) then
+         zcondicea_2(ig,l)= -zfallice_2(ig,l+1)
+       end if
+       zfallice_2(ig,l) = zcondicea_2(ig,l)+zfallice_2(ig,l+1)
+     end if
+   end do
+ end do
+ diff_zcondicea(:,:) = zcondicea_2(:,:) - zcondicea(:,:)
+ diff_zfallice(:) = zfallice_2(:,1) - zfallice(:)
+ call writediagfi(ngrid, "diff_zfallice", "sans - avec microphysique", "", 2, diff_zfallice)
+ call writediagfi(ngrid, "diff_zcondicea", "sans - avec microphysique", "", 3, diff_zcondicea)
 !----------------------------------------------------------------------------------------------------------------------!
 ! 5. Main co2condens
@@ -307,9 +364,5 @@
 !   5.3. Compute zfallheat
 !----------------------------------------------------------------------------------------------------------------------!
-      if (zfallice(ig)>0) then
-        zfallheat = zfallice(ig) * ( pphi(ig,1) - phisfi(ig) + cpice*(ztcond(ig,1) - ztcondsol(ig)) )/latcond
-      else
-        zfallheat = 0.
-      end if
+      zfallheat = 0.
 !----------------------------------------------------------------------------------------------------------------------!
 !   5.4. Compute direct condensation/sublimation of CO2 ice
@@ -320,12 +373,10 @@
 !   5.4.a. If there is not enough CO2 tracer in 1st layer to condense
 !----------------------------------------------------------------------------------------------------------------------!
-      if (igcm_co2/=0) then
-!       Available CO2 tracer in layer 1 at end of timestep (kg/m2)
-        availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) + (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)*ptimestep) )
-        if ( zcondices(ig) * ptimestep>availco2 ) then
-          zcondices(ig) = availco2/ptimestep
-          zdiceco2(ig) = zcondices(ig) + zfallice(ig)
-          pdtsrfc(ig) = (latcond/pcapcal(ig)) *  (zcondices(ig)+zfallheat)
-        end if
+!     Available CO2 tracer in layer 1 at end of timestep (kg/m2)
+      availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) + (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)*ptimestep) )
+      if ( zcondices(ig) * ptimestep>availco2 ) then
+        zcondices(ig) = availco2/ptimestep
+        zdiceco2(ig) = zcondices(ig) + zfallice(ig)
+        pdtsrfc(ig) = (latcond/pcapcal(ig)) *  (zcondices(ig)+zfallheat)
       end if
 !----------------------------------------------------------------------------------------------------------------------!
@@ -364,8 +415,4 @@
         end if
       end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.6.b. Set albedo and emissivity of the surface
-!----------------------------------------------------------------------------------------------------------------------!
-      call albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
 !----------------------------------------------------------------------------------------------------------------------!
 !   5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow)
@@ -481,11 +528,6 @@
 #else
       do l = 1, nlayer
-        if (.not. co2clouds) then
-!         Tendencies on T
-          pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l))  - zmflux(l+1)*(ztm(l+1) - zt(l)) &
-                         + zcondicea(ig,l)*(ztcond(ig,l)-zt(l)) )
-        else
-          pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l))  - zmflux(l+1)*(ztm(l+1) - zt(l)) )
-        end if
+!       Tendencies on T
+        pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l))  - zmflux(l+1)*(ztm(l+1) - zt(l)) )
 
 !       Tendencies on U
@@ -513,7 +555,11 @@
   end do  ! loop on ig
 !----------------------------------------------------------------------------------------------------------------------!
+!   5.6.b. Set albedo and emissivity of the surface
+!----------------------------------------------------------------------------------------------------------------------!
+  call albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
+!----------------------------------------------------------------------------------------------------------------------!
 ! 6. CO2 snow / clouds scheme
 !----------------------------------------------------------------------------------------------------------------------!
-  call co2snow(ngrid, nlayer, ptimestep, emisref, condsub, pplev, zcondicea, zcondices, zdqssed_co2*ptimestep, &
+  call co2snow(ngrid, nlayer, ptimestep, emisref, condsub, pplev, zcondicea, zcondices, zdqssed_co2, &
                pemisurf)
 !----------------------------------------------------------------------------------------------------------------------!
Index: trunk/LMDZ.MARS/libf/phymars/co2snow.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 2447)
@@ -31,5 +31,4 @@
 #include "callkeys.h"
 
-
 c=======================================================================
 c     Variables
@@ -50,5 +49,5 @@
                                   !   (kg/m2/s)        
      &   pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
-                           
+
       logical, intent(in) ::
      &   condsub(ngrid) ! true if there is CO2 condensation or
@@ -107,4 +106,5 @@
       do ig = 1, ngrid
         if (condsub(ig)) then
+
           if (latitude(ig).lt.0.) then
               icap = 2 ! Southern hemisphere
@@ -123,5 +123,4 @@
      &       - pemisurf(ig)
      &      ) / ptimestep
-          
           pemisurf(ig) = pemisurf(ig) + zdemisurf * ptimestep 
           if (pemisurf(ig).lt.0.1) then
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2447)
@@ -565,4 +565,10 @@
          write(*,*) " hdofrac = ",hdofrac
 
+! Activeco2ice
+         write(*,*) "Radiatively active transported atmospheric ",
+     &              "Co2 ice ?"
+         activeco2ice=.false. ! default value
+         call getin_p("activeco2ice",activeco2ice)
+         write(*,*) " activeco2ice = ",activeco2ice
 
 ! sub-grid cloud fraction: fixed
@@ -650,4 +656,10 @@
 
          
+         if (activeco2ice.and..not.co2clouds) then
+          print*,'if activeco2ice is used, co2clouds should be used too'
+          call abort_physic(modname,
+     &          "activeco2ice requires co2clouds",1)
+         endif
+
 ! water ice clouds effective variance distribution for sedimentaion       
         write(*,*) "Sed effective variance for water ice clouds ?"
@@ -834,28 +846,56 @@
          name_iaer(1) = "dust_conrath"   !! default choice is good old Conrath profile
          IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
+
          if (nq.gt.1) then
-          ! trick to avoid problems compiling with 1 tracer
-          ! and picky compilers who know name_iaer(2) is out of bounds
-          j=2
-         IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) name_iaer(2) = 
-     &                       "stormdust_doubleq" !! storm dust two-moment scheme
-         IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind)
-     &              name_iaer(3) = "stormdust_doubleq" 
-         IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) name_iaer(2) = 
-     &                       "topdust_doubleq" !! storm dust two-moment scheme
-         IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm)
-     &              name_iaer(3) =  "topdust_doubleq"
-         IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN 
-             name_iaer(2) = "stormdust_doubleq"
-             name_iaer(3) = "topdust_doubleq"
-         ENDIF
-         IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN 
-             name_iaer(3) = "stormdust_doubleq"
-             name_iaer(4) = "topdust_doubleq"
-         ENDIF
-         IF (water.AND.activice) name_iaer(j) = "h2o_ice"      !! radiatively-active clouds
-         IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
+           ! trick to avoid problems compiling with 1 tracer
+           ! and picky compilers who know name_iaer(2) is out of bounds
+           j=2
+           IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) then
+             name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme
+             j = j+1
+           END IF
+
+           IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind) then
+             name_iaer(j) = "stormdust_doubleq" 
+             j = j+1
+           END IF
+
+           IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) then
+             name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme
+             j = j+1
+           END IF
+ 
+           IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm) then
+             name_iaer(j) =  "topdust_doubleq"
+             j = j+1
+           END IF
+
+           IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN 
+             name_iaer(j) = "stormdust_doubleq"
+             name_iaer(j+1) = "topdust_doubleq"
+             j = j+2
+           ENDIF
+
+           IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN 
+             name_iaer(j) = "stormdust_doubleq"
+             name_iaer(j+1) = "topdust_doubleq"
+             j = j+2
+           ENDIF
+
+           IF (water.AND.activice) then
+            name_iaer(j) = "h2o_ice"      !! radiatively-active clouds
+            j = j+1
+           END IF
+
+           IF (co2clouds.AND.activeco2ice) then
+             name_iaer(j) = "co2_ice" !! radiatively-active co2 clouds
+             j = j+1
+           ENDIF
+
+           IF (submicron.AND.active) then
+             name_iaer(j) = "dust_submicron" !! JBM experimental stuff
+             j = j+1
+           ENDIF
          endif ! of if (nq.gt.1)
-
 c        ----------------------------------------------------------
 
Index: trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90	(revision 2447)
@@ -291,5 +291,4 @@
   ! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
   zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep
-
   where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold
 
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2447)
@@ -251,4 +251,8 @@
       real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
       real rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
+      real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius
+      real rhocloudco2(ngrid,nlayer)  ! CO2 Cloud density (kg.m-3)
+      real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the 
+                                  !   size distribution
       REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
                                       ! thermal inertia (J.s-1/2.m-2.K-1)
@@ -506,7 +510,5 @@
 
 c=======================================================================
-c---------- begin added by CM
       pdq(:,:,:) = 0.
-c---------end added by CM
 
 
@@ -544,5 +546,4 @@
      &         mem_Mccn_co2,mem_Nccn_co2,
      &         mem_Mh2o_co2,watercap)
-
 
 #else
@@ -864,5 +865,5 @@
            IF(callnlte) then
               if(nltemodel.eq.0.or.nltemodel.eq.1) then
-                 CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
+                CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
               else if(nltemodel.eq.2) then
                 co2vmr_gcm(1:ngrid,1:nlayer)=
@@ -934,8 +935,7 @@
      &     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
      &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
-     &     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
-     &     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
-     &     clearsky,totcloudfrac)
-
+     &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
+     &     rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
+     &     alpha_hmons,nohmons,clearsky,totcloudfrac)
            ! case of sub-grid water ice clouds: callradite for the clear case
             IF (CLFvarying) THEN
@@ -944,5 +944,5 @@
                ! if
                ! CLFvarying ...) ?? AP ??
-               clearsky=.true. !
+               clearsky=.true.  
                CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
      &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
@@ -952,5 +952,6 @@
      &              tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
      &              taucloudtesclf,rdust,
-     &              rice,nuice,co2ice,rstormdust,rtopdust,totstormfract,
+     &              rice,nuice,riceco2, nuiceco2,co2ice,rstormdust,
+     &              rtopdust,totstormfract,
      &              clearatm,dsords,dsotop,alpha_hmons,nohmons,
      &              clearsky,totcloudfrac)
@@ -1318,5 +1319,4 @@
      &        zcondicea_co2microp,sensibFlux,
      &        dustliftday,local_time,watercap,dwatercap_dif)
-
           DO ig=1,ngrid
              zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
@@ -1625,5 +1625,5 @@
      &           pq,pdq,zdqcloudco2,zdtcloudco2,
      &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
-     &           rsedcloud,rhocloud,zzlev,zdqssed_co2,
+     &           rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
      &           pdu,pu,zcondicea_co2microp, co2ice)
 
@@ -2156,22 +2156,4 @@
        ENDDO
 
-      ! test for co2 conservation with co2 microphysics
-      if (igcm_co2_ice.ne.0) then
-        ! Total mass of co2
-        mtotco2 = 0.
-        icetotco2 = 0.
-        vaptotco2 = 0.
-        do ig=1,ngrid
-          do l=1,nlayer
-            vaptotco2 = vaptotco2 + zq(ig,l,igcm_co2) *
-     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
-            icetotco2 = icetotco2 + zq(ig,l,igcm_co2_ice) *
-     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
-          end do
-          mtotco2 = mtotco2 + co2ice(ig)
-        end do
-        mtotco2 = mtotco2 + vaptotco2 + icetotco2
-      endif ! of if (igcm_co2_ice.ne.0)
-
 c    Compute surface stress : (NB: z0 is a common in surfdat.h)
 c     DO ig=1,ngrid
@@ -2432,5 +2414,4 @@
 c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
              enddo
-
              call watersat(ngrid*nlayer,zt,zplay,zqsat)
              satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
@@ -2818,6 +2799,6 @@
 c        WRITEDIAGFI can ALSO be called from any other subroutines
 c        for any variables !!
-c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
-c    &                  emis)
+c         call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+c     &                  emis)
 c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
 c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
@@ -2868,5 +2849,4 @@
      &                         '',2,taucloudtes)
              endif
-
 #else
       !!! this is to ensure correct initialisation of mesoscale model
@@ -3774,4 +3754,7 @@
      &                   'temperature variation of CO2 latent heat',
      &                   'K/s',3,zdtcloudco2)
+         call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     &                  emis)
+
       end if
          
@@ -3794,5 +3777,4 @@
 
       END IF       ! if(ngrid.ne.1)
-
 ! test for co2 conservation with co2 microphysics
       if (igcm_co2_ice.ne.0) then
@@ -3810,6 +3792,4 @@
      &                     'Total CO2 mass conservation in physic',
      &                     '%',0,(co2totA-co2totB)/co2totA)
-        call test_vmr(pq(:,:,igcm_co2), pdq(:,:,igcm_co2), ptimestep,
-     &                 ngrid, nlayer, 'end physiq_mod')
       endif ! of if (igcm_co2_ice.ne.0)
 ! XIOS outputs
Index: trunk/LMDZ.MARS/libf/phymars/suaer.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/suaer.F90	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/suaer.F90	(revision 2447)
@@ -154,4 +154,15 @@
 !       Reference wavelength in the infrared:
         longrefir(iaer)=12.1E-6  ! 825cm-1 TES/MGS
+!==================================================================
+        CASE("co2_ice") aerkind             ! CO2 ice crystals
+!==================================================================
+!       Visible domain:
+        file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat' !
+!       Infrared domain:
+        file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'  !
+!       Reference wavelength in the visible:
+        longrefvis(iaer) = 0.67E-6  ! 1.5um OMEGA/MEx
+!       Reference wavelength in the infrared:
+        longrefir(iaer) = 1E-6   ! 1 micron a vérifier si on prends cette valeur
 !==================================================================
         CASE("stormdust_doubleq") aerkind   ! Two-moment scheme for stormdust - radiative properties
Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 2443)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 2447)
@@ -7,10 +7,12 @@
       SUBROUTINE updatereffrad(ngrid,nlayer,
      &                rdust,rstormdust,rtopdust,rice,nuice,
-     &                reffrad,nueffrad,
+     &                reffrad,nueffrad, riceco2, nuiceco2,
      &                pq,tauscaling,tau,pplay)
        USE updaterad, ONLY: updaterdust, updaterice_micro,
-     &                      updaterice_typ
+     &                      updaterice_microco2, updaterice_typ
        use tracer_mod, only: nqmx, igcm_dust_mass, igcm_dust_number,
      &                       igcm_h2o_ice, igcm_ccn_mass, radius,
+     &                       igcm_co2_ice, nuiceco2_ref,
+     &                       igcm_ccnco2_number, igcm_ccnco2_mass,
      &                       igcm_ccn_number, nuice_ref, varian,
      &                       ref_r0, igcm_dust_submicron,
@@ -68,5 +70,9 @@
       REAL, INTENT(out) :: nueffrad(ngrid,nlayer,naerkind)
       REAL, INTENT(in) :: tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
-      
+c     CO2 ice mean radius (m)
+      REAL, INTENT(out) :: riceco2(ngrid,nlayer) ! co2 ice radius
+      REAL, INTENT(out) :: nuiceco2(ngrid,nlayer)
+
+     
 c     Local variables:
 c     ---------------
@@ -85,4 +91,6 @@
       REAL Mo,No                       ! Mass and number of ccn
       REAL rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
+c     For CO2 microphysics only:
+      REAL :: rhocloudco2(ngrid, nlayer) ! co2 cloud density
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -176,4 +184,31 @@
        ENDIF ! of if (water.AND.activice)
 
+
+c       1.3 CO2-ice particles --- in progress.....
+c       ----------------------- j'ai juste copié-collé la partie 1.2
+
+        IF (co2clouds.AND.activeco2ice) THEN
+
+c    At firstcall, the true number and true mass of cloud condensation nuclei are not known.
+c    Indeed it is scaled on the prescribed dust opacity via a 'tauscaling' coefficient
+c    computed after radiative transfer. If tauscaling is not in startfi, we make an assumption for its value.
+
+          IF (firstcall) THEN
+            !IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
+            !IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
+            firstcall = .false.
+          ENDIF
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              call updaterice_microco2(pq(ig,l,igcm_co2_ice),
+     &                              pq(ig,l,igcm_ccnco2_mass),
+     &                              pq(ig,l,igcm_ccnco2_number),
+     &                              tauscaling(ig),riceco2(ig,l),
+     &                              rhocloudco2(ig,l))
+              nuiceco2(ig,l) = nuiceco2_ref
+            END DO
+          ENDDO
+        ENDIF ! of if (co2clouds.AND.activeco2ice)
+
 c==================================================================
 c 2. Radius used in the radiative transfer code (reffrad)
@@ -227,4 +262,13 @@
           ENDDO
 c==================================================================
+        CASE("co2_ice") aerkind             ! CO2 ice crystals
+c==================================================================
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              reffrad(ig,l,iaer)=riceco2(ig,l)*(1.+nuiceco2_ref)
+              nueffrad(ig,l,iaer)=nuiceco2_ref
+            ENDDO
+          ENDDO
+c==================================================================
         CASE("stormdust_doubleq") aerkind! Two-moment scheme for
 c       stormdust; same distribution than normal dust
