Index: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2562)
@@ -159,4 +159,7 @@
       INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
       INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
+      INTEGER,SAVE :: iccnco2_h2o_number ! index of tracer containing CCN number
+      INTEGER,SAVE :: iccnco2_h2o_mass_ice ! index of tracer containing CCN number
+      INTEGER,SAVE :: iccnco2_h2o_mass_ccn ! index of tracer containing CCN number
       INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
 
@@ -241,4 +244,7 @@
           iccnco2_mass=0
           iccnco2_number=0
+          iccnco2_h2o_mass_ice=0
+          iccnco2_h2o_mass_ccn=0
+          iccnco2_h2o_number=0
           ico2_ice=0
           do iq=1,nq
@@ -254,4 +260,19 @@
               iccnco2_number=iq
               write(*,*)"callsedim: iccnco2_number=",iccnco2_number
+            endif
+            if (noms(iq).eq."ccnco2_h2o_number") then
+              iccnco2_h2o_number=iq
+              write(*,*)"callsedim: iccnco2_h2o_number=",
+     &                   iccnco2_h2o_number
+            endif
+            if (noms(iq).eq."ccnco2_h2o_mass_ice") then
+              iccnco2_h2o_mass_ice=iq
+              write(*,*)"callsedim: iccnco2_h2o_mass_ice=",
+     &                 iccnco2_h2o_mass_ice
+            endif
+            if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
+              iccnco2_h2o_mass_ccn=iq
+              write(*,*)"callsedim: iccnco2_h2o_mass_ccn=",
+     &                 iccnco2_h2o_mass_ccn
             endif
           enddo
@@ -263,4 +284,14 @@
             call abort_physic(modname,"missing co2 ccn tracers",1)
           endif
+          if (co2useh2o) then
+            if((iccnco2_h2o_mass_ccn.eq.0).or.
+     &         (iccnco2_h2o_mass_ice.eq.0).or.
+     &         (iccnco2_h2o_number.eq.0)) then
+            write(*,*) 'callsedim: error! could not identify'
+            write(*,*) ' tracers for ccn co2 mass and number mixing'
+            write(*,*) ' ratio and co2clouds are activated!'
+            call abort_physic(modname,"missing co2 ccn tracers",1)
+           end if
+          end if
        ENDIF                    !of if (co2clouds)
 
@@ -398,5 +429,8 @@
         if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
      &        (iq .ne. iccnco2_mass) .and. (iq .ne. 
-     &        iccnco2_number) .and. ! no sedim for gaz or CO2 clouds  (done in microtimestep)
+     &        iccnco2_number) .and.
+     &        (iq.ne.iccnco2_h2o_number).and.
+     &        (iq.ne.iccnco2_h2o_mass_ice).and.
+     &        (iq.ne.iccnco2_h2o_mass_ccn).and.  ! no sedim for gaz or CO2 clouds  (done in microtimestep)
      &        iq .ne. igcm_hdo_ice) then !MVals: hdo is transported by h2o
 c -----------------------------------------------------------------
Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F90	(revision 2562)
@@ -7,17 +7,8 @@
 ! Contains subroutines:
 !      - co2cloud: of co2 cloud microphysics
-!
-!      - ini_co2cloud: (only used in phys_state_var_init_mod.F90)
-!
-!      - end_co2cloud: (only used in phys_state_var_init_mod.F90)
 !======================================================================================================================!
 module co2cloud_mod
 
 implicit none
-
-double precision, allocatable, save :: &
-   mem_Mccn_co2(:,:), &! Memory of CCN mass of H2O and dust used by CO2
-   mem_Mh2o_co2(:,:), &! Memory of H2O mass integred into CO2 crystal
-   mem_Nccn_co2(:,:)   ! Memory of CCN number of H2O and dust used by CO2
 
 contains
@@ -98,6 +89,6 @@
 !======================================================================================================================!
   subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
-                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev, pdqs_sedco2, pdu, &
-                      pu, pcondicea, co2ice)
+                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev,&
+                      pdqs_sedco2, pdqs_sedccn, pdu, pu, pcondicea, co2ice)
 
   use ioipsl_getincom, only: getin
@@ -107,5 +98,6 @@
   use conc_mod,        only: mmean, rnew
   use tracer_mod,      only: igcm_co2, igcm_co2_ice, igcm_dust_mass, igcm_dust_number, igcm_h2o_ice, &
-                             igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, rho_dust, &
+                             igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, &
+                             igcm_ccnco2_h2o_number, igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, rho_dust, &
                              nuiceco2_sed, nuiceco2_ref, r3n_q, rho_ice, nuice_sed
 
@@ -159,4 +151,5 @@
      rhocloud(ngrid,nlay),      & ! Water Cloud density (kg.m-3)
      pdqs_sedco2(ngrid),        & ! CO2 flux at the surface
+     pdqs_sedccn(ngrid,nq),      & ! CCN flux at the surface
      pdqcloudco2(ngrid,nlay,nq),& ! Tendency due to CO2 condensation (kg/kg.s-1)
      pcondicea(ngrid,nlay),     & ! Rate of condensation/sublimation of co2 ice in layers
@@ -240,5 +233,6 @@
      subpdqsed(ngrid,nlay,nq),      &! Tendancy due to sedimentation
      sum_subpdqs_sedco2(ngrid),     &! CO2 flux at the surface
-!  ---For sub grid T distribution
+     sum_subpdqs_sedccn(ngrid,nq),  &!
+      !  ---For sub grid T distribution
      zt(ngrid,nlay),                &! Local value of temperature
      zq_co2vap(ngrid, nlay),        &! Local value of CO2 vap
@@ -251,4 +245,6 @@
      Nccnco2,                       &! buffer: number of ccn used for co2 condensation
      Qccnco2,                       &! buffer: mass of ccn used for co2 condensation
+     Nccnco2_h2o, &
+     Qccnco2_h2o, &
      Niceco2,                       &! buffer: mmr co2 ice
      epaisseur(ngrid,nlay),         &! Layer thickness (m)
@@ -420,4 +416,5 @@
   zqsed0(1:ngrid,1:nlay,1:nq) = 0.
   sum_subpdqs_sedco2(1:ngrid) = 0.
+  sum_subpdqs_sedccn(1:ngrid,1:nq) = 0.
   subpdqsed(1:ngrid,1:nlay,1:nq) = 0.
 !----------------------------------------------------------------------------------------------------------------------!
@@ -600,4 +597,13 @@
 
          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_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
+                                                  pdq(ig,l,igcm_ccnco2_h2o_number)
+
+        sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
+                                                    pdq(ig,l,igcm_ccnco2_h2o_mass_ice)
+
+        sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
+                                                    pdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
        end if
      end do ! ngrid
@@ -607,24 +613,24 @@
 !----------------------------------------------------------------------------------------------------------------------!
    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)
+                          subpdqcloudco2, subpdtcloudco2, nq, tauscaling, 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
+         <= 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
+       if(pq(ig,l,igcm_co2) + microtimestep*(sum_subpdq(ig,l,igcm_co2)+subpdqcloudco2(ig,l,igcm_co2)) <= 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. &
+          *microtimestep)<=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
+          *microtimestep<=1e-20)) then
          subpdqcloudco2(ig,l,igcm_ccnco2_number) = - pq(ig,l,igcm_ccnco2_number)/microtimestep + 1. &
                                                   - sum_subpdq(ig,l,igcm_ccnco2_number)
@@ -635,4 +641,29 @@
          subpdqcloudco2(ig,l,igcm_dust_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_mass)
        end if
+
+   ! ccnco2_h2o_number and masses
+       if (co2useh2o) then
+         if (((pq(ig,l,igcm_ccnco2_h2o_number) + (sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
+           subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number)) &
+            *microtimestep)<=1.) .or. &
+           (pq(ig,l,igcm_ccnco2_h2o_mass_ice)+pq(ig,l,igcm_ccnco2_h2o_mass_ccn) +&
+           (sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) +&
+           sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) &
+            *microtimestep<=1e-20)) then
+
+           subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = - pq(ig,l,igcm_ccnco2_h2o_number)/microtimestep + 1. &
+                                                         - sum_subpdq(ig,l,igcm_ccnco2_h2o_number)
+           subpdqcloudco2(ig,l,igcm_ccn_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number)
+
+           subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = - pq(ig,l,igcm_ccnco2_h2o_mass_ice)/microtimestep + 1e-20 &
+                                                           - sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice)
+
+           subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = - pq(ig,l,igcm_ccnco2_h2o_mass_ccn)/microtimestep + 1e-20 &
+                                                           - sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
+
+           subpdqcloudco2(ig,l,igcm_ccn_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) &
+                                                - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)
+         end if
+       end if
      end do
    end do
@@ -663,4 +694,14 @@
 
           sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
+
+          sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
+                                                      subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice)
+
+          sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
+                                                      subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)
+
+          sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
+                                                    subpdqcloudco2(ig,l,igcm_ccnco2_number)
+
         end if
       end do ! ngrid
@@ -686,8 +727,17 @@
           Qccnco2 = max(zqsed(ig,l,igcm_ccnco2_mass), threshold)
 
+          Nccnco2_h2o = 0.
+          Qccnco2_h2o = 0.
+          if (co2useh2o) then
+            Nccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_number), threshold)
+            Qccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_mass_ice)+zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn), threshold)
+            Nccnco2 = Nccnco2 - Nccnco2_h2o
+            Qccnco2 = Qccnco2 - Qccnco2_h2o
+          end if
+
           ! Get density cloud and co2 ice particle radius
-          if (Niceco2.ne.0d0) then
-            call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), ztsed(ig,l), tauscaling(ig), &
-                                     riceco2(ig,l), rhocloudco2t(ig,l))
+          if (Niceco2/=0d0) then
+            call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o),&
+                                     ztsed(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
           else
             riceco2(ig,l) = 0.
@@ -702,4 +752,10 @@
       zqsed0(:,:,igcm_ccnco2_mass) = zqsed(:,:,igcm_ccnco2_mass)
       zqsed0(:,:,igcm_ccnco2_number) = zqsed(:,:,igcm_ccnco2_number)
+
+      if (co2useh2o) then
+        zqsed0(:,:,igcm_ccnco2_h2o_number) = zqsed(:,:,igcm_ccnco2_h2o_number)
+        zqsed0(:,:,igcm_ccnco2_h2o_mass_ice) = zqsed(:,:,igcm_ccnco2_h2o_mass_ice)
+        zqsed0(:,:,igcm_ccnco2_h2o_mass_ccn) = zqsed(:,:,igcm_ccnco2_h2o_mass_ccn)
+      end if
 !----------------------------------------------------------------------------------------------------------------------!
 ! 4.3.c. Sedimentation of co2 ice
@@ -726,5 +782,7 @@
       call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
                    rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_mass), wq, beta)
-      !TODO: ajouter le calcule de la tendance a la surface comme co2ice
+      do ig = 1, ngrid
+        sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) = sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) + wq(ig,1) / microtimestep
+      end do
 
       wq(:,:) = 0.
@@ -732,5 +790,36 @@
       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
+      do ig = 1, ngrid
+        sum_subpdqs_sedccn(ig,igcm_ccnco2_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_number) + wq(ig,1) / microtimestep
+      end do
+
+      ! for ccnco2_h2o_mass_ice
+      if (co2useh2o) then
+        wq(:,:) = 0.
+        call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
+                     rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ice), wq, beta)
+        do ig = 1, ngrid
+          sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) + &
+                                                            wq(ig,1) / microtimestep
+        end do
+
+        wq(:,:) = 0.
+        ! for ccnco2_h2o_mass_ccn
+        call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
+                     rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ccn), wq, beta)
+        do ig = 1, ngrid
+          sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) + &
+                                                            wq(ig,1) / microtimestep
+        end do
+
+        wq(:,:) = 0.
+        ! for ccnco2_h2o_number
+        call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
+                     rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_number), wq, beta)
+        do ig = 1, ngrid
+          sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) + &
+                                                            wq(ig,1) / microtimestep
+        end do
+      end if
 !----------------------------------------------------------------------------------------------------------------------!
 ! 4.3.e. Compute tendencies due to the sedimation process
@@ -745,4 +834,15 @@
 
           subpdqsed(ig,l,igcm_co2_ice) = ( zqsed(ig,l,igcm_co2_ice) - zqsed0(ig,l,igcm_co2_ice) ) / microtimestep
+
+          if (co2useh2o) then
+          subpdqsed(ig,l,igcm_ccnco2_h2o_number) = ( zqsed(ig,l,igcm_ccnco2_h2o_number) - &
+                                                     zqsed0(ig,l,igcm_ccnco2_h2o_number) ) / microtimestep
+
+          subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ice) - &
+                                                     zqsed0(ig,l,igcm_ccnco2_h2o_mass_ice) ) / microtimestep
+
+          subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
+                                                     zqsed0(ig,l,igcm_ccnco2_h2o_mass_ccn) ) / microtimestep
+          end if
         end do
       end do
@@ -755,4 +855,14 @@
 
           sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqsed(ig,l,igcm_co2_ice)
+          if (co2useh2o) then
+            sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
+                                                        subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice)
+
+            sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
+                                                        subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn)
+
+            sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
+                                                        subpdqsed(ig,l,igcm_ccnco2_h2o_number)
+          end if
         end do
       end do
@@ -772,5 +882,7 @@
   do ig = 1, ngrid
     pdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) / real(imicroco2)
+    pdqs_sedccn(ig,:) = sum_subpdqs_sedccn(ig,:) / real(imicroco2)
   end do
+
   ! temperature tendency (T.s-1)
   do l = 1, nlay
@@ -808,4 +920,14 @@
         pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
                                             pdq(ig,l,igcm_ccn_number)
+
+        pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_number) / real(imicroco2) ) - &
+                                                   pdq(ig,l,igcm_ccnco2_h2o_number)
+
+        pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) / real(imicroco2) ) - &
+                                                   pdq(ig,l,igcm_ccnco2_h2o_mass_ice)
+
+        pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) / real(imicroco2) )- &
+                                                        pdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
+
       end if
     end do ! ngrid
@@ -831,10 +953,25 @@
       myT = pteff(ig,l) + (pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
 
+      Nccnco2_h2o = 0.
+      Qccnco2_h2o = 0.
+      if (co2useh2o) then
+      Nccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_number) + (pdq(ig,l,igcm_ccnco2_h2o_number) + &
+                       pdqcloudco2(ig,l, igcm_ccnco2_h2o_number)) * ptimestep) &
+                    , threshold)
+
+      Qccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_mass_ice) + pqeff(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
+                     (pdq(ig,l,igcm_ccnco2_h2o_mass_ice) + pdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
+                      pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) + pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) * &
+                      ptimestep)&
+                    , threshold)
+      Nccnco2 = Nccnco2 - Nccnco2_h2o
+      Qccnco2 = Qccnco2 - Qccnco2_h2o
+      end if
       ! Compute particle size
-      call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), myT, tauscaling(ig), riceco2(ig,l), &
-                               rhocloudco2(ig,l))
+      call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o), myT, &
+                               tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
 
      ! Compute opacities
-      if ( (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1.) ) then
+      if ( (Niceco2 <= threshold .or. (Nccnco2)*tauscaling(ig) <= 1.) ) then ! Nccnco2 inclut Nccnco2_h2o
         riceco2(ig,l) = 0.
         Qext1bins2(ig,l) = 0.
@@ -844,8 +981,8 @@
 
         do i = 1, nbinco2_cld
-          n_aer(i) = -0.5 * Nccnco2*tauscaling(ig) * n_derf
+          n_aer(i) = -0.5 * (Nccnco2)*tauscaling(ig) * n_derf
 
           n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2)
-          n_aer(i) = n_aer(i) + (0.5 * Nccnco2*tauscaling(ig) * n_derf)
+          n_aer(i) = n_aer(i) + (0.5 * (Nccnco2)*tauscaling(ig) * n_derf)
 
           Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
@@ -928,4 +1065,10 @@
 
           pdqcloudco2(ig,l,igcm_ccn_number) = pdqcloudco2(ig,l,igcm_ccn_number) * co2cloudfrac(ig,l)
+
+          pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) * co2cloudfrac(ig,l)
+
+          pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) * co2cloudfrac(ig,l)
+
+          pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) * co2cloudfrac(ig,l)
         end if
       end do ! ngrid
@@ -960,9 +1103,4 @@
   call WRITEDIAGFI(ngrid, "tau1mic", "co2 ice opacity 1 micron", " ", 2, tau1mic)
 
-  call WRITEDIAGFI(ngrid, "mem_Nccn_co2", "CCN number used by CO2", "kg.kg-1", 3, mem_Nccn_co2)
-
-  call WRITEDIAGFI(ngrid, "mem_Mccn_co2", "CCN mass used by CO2",  "kg.kg-1", 3, mem_Mccn_co2)
-
-  call WRITEDIAGFI(ngrid, "mem_Mh2o_co2", "H2O mass in CO2 crystal", "kg.kg-1", 3, mem_Mh2o_co2)
   if (CLFvaryingCO2) then
     call WRITEDIAGFI(ngrid, "co2cloudfrac", "co2 cloud fraction", " ", 3, co2cloudfrac)
@@ -973,49 +1111,3 @@
   end subroutine co2cloud
 
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-
-!======================================================================================================================!
-! SUBROUTINE: ini_co2cloud ============================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Allocate arrays used for co2 microphysics
-!======================================================================================================================!
-  subroutine ini_co2cloud(ngrid,nlayer)
-
-  implicit none
-
-  integer, intent(in) :: &
-     ngrid, &! number of atmospheric columns
-     nlayer  ! number of atmospheric layers
-
-  allocate(mem_Nccn_co2(ngrid,nlayer))
-  allocate(mem_Mccn_co2(ngrid,nlayer))
-  allocate(mem_Mh2o_co2(ngrid,nlayer))
-
-  end subroutine ini_co2cloud
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-!======================================================================================================================!
-! SUBROUTINE: end_co2cloud ============================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Deallocate arrays used for co2 microphysics
-!======================================================================================================================!
-  subroutine end_co2cloud
-
-  implicit none
-
-  if (allocated(mem_Nccn_co2)) deallocate(mem_Nccn_co2)
-  if (allocated(mem_Mccn_co2)) deallocate(mem_Mccn_co2)
-  if (allocated(mem_Mh2o_co2)) deallocate(mem_Mh2o_co2)
-
-  end subroutine end_co2cloud
-
 end module co2cloud_mod
Index: trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90	(revision 2562)
@@ -61,6 +61,5 @@
 !======================================================================================================================!
   subroutine 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)
+                               subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2)
 
   use comcstfi_h,   only: pi, g, cpp
@@ -69,9 +68,11 @@
 
   use tracer_mod,   only: igcm_dust_mass, igcm_dust_number, rho_dust, igcm_h2o_ice, igcm_ccn_mass, igcm_ccn_number, &
-                          nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, nuiceco2_sed, &
+                          nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, &
+                          igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, igcm_ccnco2_h2o_number, nuiceco2_sed, &
                           nuiceco2_ref
 
   use conc_mod,     only: mmean
 
+  use nucleaco2_mod, only: nucleaco2
   use datafile_mod, only: datadir
 
@@ -163,4 +164,8 @@
 
   real :: &
+    Nccnco2, &
+    Qccnco2, &
+    Nccnco2_h2o, &
+    Qccnco2_h2o, &
     masse(ngrid,nlay),     &! mass layer (kg.m-2)
     rice(ngrid,nlay),      &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
@@ -189,12 +194,8 @@
     n_derf,                            &! derf( (rb_cldco2(1)+Rn) *dev3)
     m_derf,                            &! derf( (rb_cldco2(1)+Rm) *dev2)
-    mem_Mccn_co2(ngrid,nlay),          &! Memory of CCN mass of H2O and dust used by CO2
-    mem_Mh2o_co2(ngrid,nlay),          &! Memory of H2O mass integred into CO2 crystal
-    mem_Nccn_co2(ngrid,nlay),          &! Memory of CCN number of H2O and dust used by CO2
     n_aer(nbinco2_cld),                &! Radius used by the microphysical scheme (m)
     m_aer(nbinco2_cld),                &! number concentration V-1 of particle/each size bin
     n_aer_h2oice(nbinco2_cld),         &! mass mixing ratio of particle/each size bin
     m_aer_h2oice(nbinco2_cld),         &! Same - for CO2 nucleation
-    rad_h2oice(nbinco2_cld),           &! Same - for CO2 nucleation
     Ic_rice,                           &! Mass transfer rate CO2 ice crystal
     ratioh2o_ccn,                      &! 1./(zq(ig,l,igcm_h2o_ice)  + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
@@ -374,54 +375,8 @@
           end do
 
-          ! 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)**3
-             end do
-          end if
-
-          ! 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) + threshold
-
-            ! Total mass of H20 crystals,CCN included
-            No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
-
-            Rn = -log(rice(ig,l))
-
-            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)
-            end do
-          end if
-
           ! Call to nucleation routine
-          call nucleaCO2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, n_aer_h2oice, rad_h2oice, rateh2o, vo2co2)
-
+          call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, vo2co2, mtetaco2)
           dN = 0.
           dM = 0.
-          dNh2o = 0.
-          dMh2o = 0.
 
           do i = 1, nbinco2_cld
@@ -431,5 +386,60 @@
           end do
 
-          if (co2useh2o) then
+          ! Now increment CCN tracers and update dust tracers
+          dNN = min(dN,zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
+          dMM = min(dM,zq(ig,l,igcm_dust_mass))  ! idem dans le min
+
+          zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
+
+          zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
+
+          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
+
+          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
+
+          ! 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)**3
+             end do
+          end if
+        end if ! of if No > 1e-30
+
+        ! Same but with h2o particles as CCN only if co2useh2o = .true.
+        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))
+
+          ! Total mass of H20 crystals,CCN included
+          Mo = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass) * tauscaling(ig) + threshold
+
+          No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
+          if (No > threshold) then
+            Rn = -log(rice(ig,l))
+
+            Rm = Rn - 3. * sigma_ice * sigma_ice
+
+            n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
+            m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
+
+            n_aer_h2oice(:)=0.
+            m_aer_h2oice(:)=0.
+            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
+            end do
+
+            call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_h2oice, rateh2o, vo2co2, mteta)
+            dNh2o = 0.
+            dMh2o = 0.
+
             do i = 1, nbinco2_cld
               Probah2o = 1.0 - exp(-1.*microtimestep*rateh2o(i))
@@ -437,20 +447,6 @@
               dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
             end do
-          end if
-
-          ! Now increment CCN tracers and update dust tracers
-          dNN = min(dN,zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
-          dMM = min(dM,zq(ig,l,igcm_dust_mass))  ! idem dans le min
-
-          zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
-
-          zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
-
-          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
-
-          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
-
-          ! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it
-          if (co2useh2o) then
+
+            ! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it
             dNNh2o = dNh2o/tauscaling(ig)
             dNNh2o = min(dNNh2o,zq(ig,l,igcm_ccn_number))
@@ -465,19 +461,16 @@
             dMh2o_ice = min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
 
-            zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice + dMh2o_ccn
-
-            zq(ig,l,igcm_ccnco2_number) =  zq(ig,l,igcm_ccnco2_number) + dNNh2o
-
-            zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number)  - dNNh2o
-
+            zq(ig,l,igcm_ccnco2_h2o_mass_ice) = zq(ig,l,igcm_ccnco2_h2o_mass_ice)  + dMh2o_ice
+            zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = zq(ig,l,igcm_ccnco2_h2o_mass_ccn)  + dMh2o_ccn
+            zq(ig,l,igcm_ccnco2_h2o_number) = zq(ig,l,igcm_ccnco2_h2o_number) + dNNh2o
+
+            zq(ig,l,igcm_ccn_number) =  zq(ig,l,igcm_ccn_number)  - dNNh2o
             zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)  - dMh2o_ice
-
             zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dMh2o_ccn
 
-            mem_Mh2o_co2(ig,l) = mem_Mh2o_co2(ig,l) + dMh2o_ice
-            mem_Mccn_co2(ig,l) = mem_Mccn_co2(ig,l) + dMh2o_ccn
-            mem_Nccn_co2(ig,l) = mem_Nccn_co2(ig,l) + dNNh2o
-          end if ! of if co2useh2o
-        end if ! of if No > 1e-30
+	          zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice + dMh2o_ccn
+	          zq(ig,l,igcm_ccnco2_number) =  zq(ig,l,igcm_ccnco2_number) + dNNh2o
+          end if
+        end if ! of if co2useh2o
       end if ! of is satu > 1
 !----------------------------------------------------------------------------------------------------------------------!
@@ -489,9 +482,19 @@
 !----------------------------------------------------------------------------------------------------------------------!
       ! we trigger crystal growth
-      if (zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) + threshold >= 1) then
-
-        call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(zq(ig,l,igcm_ccnco2_mass)), &
-                                 dble(zq(ig,l,igcm_ccnco2_number)), zt(ig,l), tauscaling(ig), riceco2(ig,l), &
-                                 rhocloudco2(ig,l))
+      if ((zq(ig,l,igcm_ccnco2_number)) * tauscaling(ig) + threshold >= 1) then
+        Nccnco2 = dble(zq(ig,l,igcm_ccnco2_number))
+        Qccnco2 = dble(zq(ig,l,igcm_ccnco2_mass))
+        Nccnco2_h2o = 0.
+        Qccnco2_h2o = 0.
+
+        if (co2useh2o) then
+          Nccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_number)
+          Qccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
+          Nccnco2 = Nccnco2 - Nccnco2_h2o
+          Qccnco2 = Qccnco2 - Qccnco2_h2o
+        end if
+
+        call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), &
+                                 dble(Nccnco2_h2o), zt(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
         Ic_rice = 0.
 
@@ -532,35 +535,23 @@
       ! On sublime tout
       if ((zq(ig,l,igcm_co2_ice) <= threshold).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
-
-        if (co2useh2o) then
-
-          if (mem_Mccn_co2(ig,l) > 0) then
-            zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + mem_Mccn_co2(ig,l)
+          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass)
+          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)
+          zq(ig,l,igcm_ccnco2_mass) = 0.
+          zq(ig,l,igcm_ccnco2_number) = 0.
+          if (co2useh2o) then
+            zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
+                                      zq(ig,l,igcm_ccnco2_h2o_mass_ice)
+            zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_h2o_number)
+
+            zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
+            zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ice)
+            zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + zq(ig,l,igcm_ccnco2_h2o_number)
+            zq(ig,l,igcm_ccnco2_h2o_number) = 0.
+            zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = 0.
+            zq(ig,l,igcm_ccnco2_h2o_mass_ice) = 0.
           end if
-
-          if (mem_Mh2o_co2(ig,l) > 0) then
-            zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + mem_Mh2o_co2(ig,l)
-          end if
-
-          if (mem_Nccn_co2(ig,l) > 0) then
-            zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + mem_Nccn_co2(ig,l)
-          end if
-
-        end if
-
-        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass) - ( mem_Mh2o_co2(ig,l) + &
-                                     mem_Mccn_co2(ig,l) )
-
-        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)  - mem_Nccn_co2(ig,l)
-
-        zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)  + zq(ig,l,igcm_co2_ice)
-
-        zq(ig,l,igcm_co2_ice) = 0.
-        zq(ig,l,igcm_ccnco2_mass) = 0.
-        zq(ig,l,igcm_ccnco2_number) = 0.
-        mem_Nccn_co2(ig,l) = 0.
-        mem_Mh2o_co2(ig,l) = 0.
-        mem_Mccn_co2(ig,l) = 0.
-        riceco2(ig,l) = 0.
+          zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) + zq(ig,l,igcm_co2_ice)
+          zq(ig,l,igcm_co2_ice) = 0.
+          riceco2(ig,l) = 0.
       end if !of if co2_ice < threshold or zq(ccnco2_number) < 1
     end do ! of ig loop
@@ -587,4 +578,13 @@
 
     subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
+
+    subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) = (zq(:,:,igcm_ccnco2_h2o_mass_ice)-zq0(:,:,igcm_ccnco2_h2o_mass_ice)&
+                                                   )/microtimestep
+
+    subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) = (zq(:,:,igcm_ccnco2_h2o_mass_ccn)-zq0(:,:,igcm_ccnco2_h2o_mass_ccn)&
+                                                   )/microtimestep
+
+    subpdqcloudco2(:,:,igcm_ccnco2_h2o_number) = ( zq(:,:,igcm_ccnco2_h2o_number) - zq0(:,:,igcm_ccnco2_h2o_number) &
+                                                 )/microtimestep
   end if
 !======================================================================================================================!
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2562)
@@ -64,4 +64,7 @@
       igcm_ccnco2_mass=0
       igcm_ccnco2_number=0
+      igcm_ccnco2_h2o_mass_ice=0
+      igcm_ccnco2_h2o_mass_ccn=0
+      igcm_ccnco2_h2o_number=0
       igcm_dust_mass=0
       igcm_dust_number=0
@@ -425,4 +428,18 @@
               count=count+1
            endif
+           if (co2useh2o) then
+           if (noms(iq).eq."ccnco2_h2o_number") then
+              igcm_ccnco2_h2o_number=iq
+              count=count+1
+           endif
+           if (noms(iq).eq."ccnco2_h2o_mass_ice") then
+              igcm_ccnco2_h2o_mass_ice=iq
+              count=count+1
+           endif
+           if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
+              igcm_ccnco2_h2o_mass_ccn=iq
+              count=count+1
+           endif
+           end if
         endif
       enddo                     ! of do iq=1,nq
Index: trunk/LMDZ.MARS/libf/phymars/nucleaco2.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nucleaco2.F90	(revision 2562)
+++ trunk/LMDZ.MARS/libf/phymars/nucleaco2.F90	(revision 2562)
@@ -0,0 +1,115 @@
+module nucleaco2_mod
+  implicit none
+  contains
+subroutine nucleaco2(pco2,temp,sat,n_ccn,nucrate,vo2co2, teta)
+USE comcstfi_h
+
+implicit none
+!*                                                     *
+!*   This subroutine computes the nucleation rate      *
+!*   as given in Pruppacher & Klett (1978) in the      *
+!*   case of water ice forming on a solid substrate.   *
+!*     Definition refined by Keese (jgr,1989)          *
+!*   Authors: F. Montmessin                            *
+!*     Adapted for the LMD/GCM by J.-B. Madeleine      *
+!*     (October 2011)                                  *
+!*     Optimisation by A. Spiga (February 2012)        *
+!*   CO2 nucleation routine dev. by Constantino        *
+!*     Listowski and Joachim Audouard (2016-2017),     *
+!*     adapted from the water ice nucleation
+!* It computes two different nucleation rates : one
+!* on the dust CCN distribution and the other one on
+!* the water ice particles distribution
+!*******************************************************
+ ! nucrate = output
+ ! nucrate_h2o en sortie aussi :
+!nucleation sur dust et h2o separement ici
+
+include "microphys.h"
+include "callkeys.h"
+
+!     Inputs
+DOUBLE PRECISION, intent(in) :: pco2,sat,vo2co2
+DOUBLE PRECISION, intent(in) :: n_ccn(nbinco2_cld)
+REAL, intent(in) :: temp !temperature
+REAL, intent(in) :: teta
+!     Output
+DOUBLE PRECISION, intent(out) :: nucrate(nbinco2_cld)
+
+!     Local variables
+DOUBLE PRECISION nco2
+DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
+DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
+DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
+!DOUBLE PRECISION fshapeco2   ! function defined at the end of the file
+DOUBLE PRECISION deltaf
+double precision mtetalocal ! local teta in double precision
+double precision fshapeco2simple,zefshapeco2
+integer i
+!*************************************************
+
+mtetalocal = dble(teta)
+
+nco2   = pco2 / kbz / temp
+rstar  = 2. * sigco2 * vo2co2 / (kbz*temp*log(sat))
+gstar  = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2)
+
+fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal) / 4.
+
+!c       Loop over size bins
+do i=1,nbinco2_cld
+
+  if ( n_ccn(i) .lt. 1e-10 ) then
+    !c no dust, no need to compute nucleation!
+    nucrate(i)=0.
+  else
+    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)
+
+    if (deltaf.eq.-100.) then
+        nucrate(i) = 0.
+    else
+        nucrate(i) = dble(sqrt ( fistar / (3.*pi*kbz*temp*(gstar*gstar)) ) * kbz * temp * rstar * rstar * 4. * pi &
+            * ( nco2*rad_cldco2(i) ) * ( nco2*rad_cldco2(i) ) / ( zefshapeco2 * nusco2 * m0co2 )  * exp (deltaf))
+    endif
+  endif ! if n_ccn(i) .lt. 1e-10
+
+enddo
+
+end subroutine nucleaco2
+
+
+!*********************************************************
+double precision function fshapeco2(cost,rap)
+implicit none
+!*        function computing the f(m,x) factor           *
+!* related to energy required to form a critical embryo  *
+!*********************************************************
+
+double precision cost,rap
+double precision yeah
+
+!! PHI
+yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
+!! FSHAPECO2 = TERM A
+fshapeco2 = (1.-cost*rap) / yeah
+fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2
+fshapeco2 = 1. + fshapeco2
+!! ... + TERM B
+yeah = (rap-cost)/yeah
+fshapeco2 = fshapeco2 + rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
+!! ... + TERM C
+fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.)
+!! FACTOR 1/2
+fshapeco2 = 0.5*fshapeco2
+
+end function fshapeco2
+
+end module nucleaco2_mod
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2562)
@@ -8,6 +8,5 @@
 subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
                      day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,co2ice, &
-                     tauscaling,totcloudfrac,wstar,mem_Mccn_co2,mem_Nccn_co2,&
-                     mem_Mh2o_co2,watercap)
+                     tauscaling,totcloudfrac,wstar,watercap)
 
   use tracer_mod, only: noms ! tracer names
@@ -28,5 +27,5 @@
 !======================================================================
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
-!  Adaptation à Mars : Yann Wanherdrick 
+!  Adaptation � Mars : Yann Wanherdrick 
 ! Objet: Lecture de l etat initial pour la physique
 ! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
@@ -66,7 +65,4 @@
   real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
   real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
-  real,intent(out) :: mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
-  real,intent(out) :: mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
-  real,intent(out) :: mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
   real,intent(out) :: watercap(ngrid) ! h2o_ice_cover
 !======================================================================
@@ -363,44 +359,4 @@
             minval(co2ice), maxval(co2ice)
 
-! Memory of the origin of the co2 particles
-if (startphy_file) then
-   call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
-   if (.not.found) then
-     write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
-     mem_Mccn_co2(:,:)=0
-   endif
-else
-   mem_Mccn_co2(:,:)=0
-endif !if (startphy_file)
-write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
-write(*,*) " <mem_Mccn_co2> range:", &
-             minval(mem_Mccn_co2), maxval(mem_Mccn_co2)
-
-if (startphy_file) then
-   call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
-   if (.not.found) then
-     write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
-     mem_Nccn_co2(:,:)=0
-   endif
-else
-  mem_Nccn_co2(:,:)=0
-endif ! if (startphy_file)
-write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
-write(*,*) " <mem_Nccn_co2> range:", &
-             minval(mem_Nccn_co2), maxval(mem_Nccn_co2)
-
-if (startphy_file) then
-   call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
-   if (.not.found) then
-     write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
-     mem_Mh2o_co2(:,:)=0
-   endif
-else
-   mem_Mh2o_co2(:,:)=0
-endif ! if (startphy_file)
-write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
-write(*,*) " <mem_Mh2o_co2> range:", &
-             minval(mem_Mh2o_co2), maxval(mem_Mh2o_co2)
-
 ! Dust conversion factor
 if (startphy_file) then
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2562)
@@ -150,5 +150,5 @@
                     phystep,time,tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,&
                     tauscaling,totcloudfrac,wstar, &
-                    mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2, watercap)
+                    watercap)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -181,7 +181,4 @@
   real,intent(in) :: totcloudfrac(ngrid)
   real,intent(in) :: wstar(ngrid)
-  real,intent(in) :: mem_Mccn_co2(ngrid,nlay) ! CCN mass of H2O and dust used by CO2
-  real,intent(in) :: mem_Nccn_co2(ngrid,nlay) ! CCN number of H2O and dust used by CO2
-  real,intent(in) :: mem_Mh2o_co2(ngrid,nlay) ! H2O mass integred into CO2 crystal
   real,intent(in) :: watercap(ngrid)
   
@@ -311,11 +308,5 @@
     enddo
   endif
-  ! Memory of the origin of the co2 particles
-  if (co2useh2o) then
-     call put_field("mem_Mccn_co2","CCN mass of H2O and dust used by CO2",mem_Mccn_co2,time)
-     call put_field("mem_Nccn_co2","CCN number of H2O and dust used by CO2",mem_Nccn_co2,time)
-     call put_field("mem_Mh2o_co2","H2O mass integred into CO2 crystal",mem_Mh2o_co2,time)
-  endif
- 
+
   ! Non-orographic gavity waves
   if (calllott_nonoro) then
Index: trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2562)
@@ -48,5 +48,4 @@
       use tracer_mod, only: ini_tracer_mod,end_tracer_mod
       use time_phylmdz_mod, only: init_time
-      use co2cloud_mod, only: ini_co2cloud,end_co2cloud
       use compute_dtau_mod, only: ini_compute_dtau_mod, &
                                   end_compute_dtau_mod
@@ -129,10 +128,5 @@
       call end_turb_mod
       call ini_turb_mod(ngrid,nlayer)
-      
-      ! allocate arrays in "co2cloud" : 
-      ! Memory of the origin of the co2 particles      
-      call end_co2cloud
-      call ini_co2cloud(ngrid,nlayer)
-      
+
       ! allocate arrays in "compute_dtau_mod":
       call end_compute_dtau_mod
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2562)
@@ -19,6 +19,5 @@
       use co2condens_mod, only: co2condens
       use co2condens_mod4micro, only: co2condens4micro
-      use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
-     &                        mem_Nccn_co2
+      use co2cloud_mod, only: co2cloud
       use callradite_mod, only: callradite
       use callsedim_mod, only: callsedim
@@ -31,4 +30,7 @@
      &                      igcm_ccn_mass, igcm_ccn_number,
      &                      igcm_ccnco2_mass, igcm_ccnco2_number,
+     &                      igcm_ccnco2_h2o_mass_ice,
+     &                      igcm_ccnco2_h2o_mass_ccn,
+     &                      igcm_ccnco2_h2o_number,
      &                      rho_ice_co2,nuiceco2_sed,nuiceco2_ref,
      &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
@@ -265,4 +267,5 @@
       DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
       real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
+      real zdqssed_ccn(ngrid,nq)  ! CCN flux at the surface  (kg.m-2.s-1)
       real zcondicea_co2microp(ngrid,nlayer)
 c     Variables used by the photochemistry
@@ -551,6 +554,5 @@
      &         tsurf,tsoil,albedo,emis,
      &         q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,
-     &         mem_Mccn_co2,mem_Nccn_co2,
-     &         mem_Mh2o_co2,watercap)
+     &         watercap)
 
 #else
@@ -573,7 +575,5 @@
       !!! a couple initializations (dummy for mesoscale) done in phyetat0
       !!! --- maybe this should be done in update_inputs_physiq_mod
-      mem_Mccn_co2(:,:)=0.0
-      mem_Nccn_co2(:,:)=0.0
-      mem_Mh2o_co2(:,:)=0.0
+
       tauscaling(:)=1.0 !! probably important
       totcloudfrac(:)=1.0
@@ -1653,4 +1653,5 @@
 c               imicroco2=50     micro-timestep is 1/50 of physical timestep
         zdqssed_co2(:) = 0.
+        zdqssed_ccn(:,:) = 0.
 
          IF (co2clouds) THEN
@@ -1660,11 +1661,16 @@
      &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
      &           rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
-     &           pdu,pu,zcondicea_co2microp, co2ice)
-
+     &           zdqssed_ccn,pdu,pu,zcondicea_co2microp, co2ice)
+
+          DO iq=1, nq
+            DO ig=1,ngrid
+              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqssed_ccn(ig,iq)
+            ENDDO  ! (ig)
+           ENDDO    ! (iq)
 c Temperature variation due to latent heat release
-               pdt(1:ngrid,1:nlayer) = 
+               pdt(1:ngrid,1:nlayer) =
      &              pdt(1:ngrid,1:nlayer) + 
      &              zdtcloudco2(1:ngrid,1:nlayer)
-            
+
 
 ! increment dust and ccn masses and numbers
@@ -1701,4 +1707,16 @@
      &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
      &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
+
+               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
+     &                pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
+     &                zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice)
+
+               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
+     &                pdq(:,:,igcm_ccnco2_h2o_mass_ccn) +
+     &                zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn)
+
+               pdq(:,:,igcm_ccnco2_h2o_number) =
+     &                   pdq(:,:,igcm_ccnco2_h2o_number) +
+     &                   zdqcloudco2(:,:,igcm_ccnco2_h2o_number)
 
 c     Negative values?
@@ -1718,4 +1736,33 @@
      &              - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
                 end where
+             where (pq(:,:,igcm_ccnco2_h2o_mass_ice) +
+     &              pq(:,:,igcm_ccnco2_h2o_mass_ccn) +
+     &              (pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
+     &               pdq(:,:,igcm_ccnco2_h2o_mass_ccn)
+     &               )*ptimestep < 0.)
+               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
+     &            - pq(:,:,igcm_ccnco2_h2o_mass_ice)
+     &               /ptimestep + 1.e-30
+               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
+     &            - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
+     &               /ptimestep + 1.e-30
+               pdq(:,:,igcm_ccnco2_h2o_number) =
+     &            - pq(:,:,igcm_ccnco2_h2o_number)
+     &                /ptimestep + 1.e-30
+             end where
+
+             where (pq(:,:,igcm_ccnco2_h2o_number) +
+     &              (pdq(:,:,igcm_ccnco2_h2o_number)
+     &               )*ptimestep < 0.)
+               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
+     &            - pq(:,:,igcm_ccnco2_h2o_mass_ice)
+     &               /ptimestep + 1.e-30
+               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
+     &            - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
+     &               /ptimestep + 1.e-30
+               pdq(:,:,igcm_ccnco2_h2o_number) =
+     &            - pq(:,:,igcm_ccnco2_h2o_number)
+     &                /ptimestep + 1.e-30
+             end where
              endif ! of if (co2useh2o)
 c     Negative values?
@@ -1960,5 +2007,6 @@
      $                  zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
      $                  fluxsurf_sw,zls,
-     $                  zdqssed_co2,zcondicea_co2microp)
+     $                  zdqssed_co2,zcondicea_co2microp
+     &                  )
          ! no scavenging yet
          zdqsc(:,:) = 0.
@@ -2366,5 +2414,5 @@
      .                tsurf,tsoil,co2ice,albedo,emis,
      .                q2,qsurf,tauscaling,totcloudfrac,wstar,
-     .                mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2,watercap)
+     .                watercap)
           
          ENDIF
Index: trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2562)
@@ -49,4 +49,7 @@
       integer,save :: igcm_ccnco2_mass   ! CCN (dust and/or water ice) for CO2 mass mixing ratio
       integer,save :: igcm_ccnco2_number ! CCN (dust and/or water ice) for CO2 number mixing ratio
+      integer,save :: igcm_ccnco2_h2o_mass_ice   ! CCN (dust and/or water ice) for CO2 mass mixing ratio
+      integer,save :: igcm_ccnco2_h2o_mass_ccn   ! CCN (dust and/or water ice) for CO2 mass mixing ratio
+      integer,save :: igcm_ccnco2_h2o_number ! CCN (dust and/or water ice) for CO2 number mixing ratio
 
       ! water
Index: trunk/LMDZ.MARS/libf/phymars/updaterad.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updaterad.F90	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/updaterad.F90	(revision 2562)
@@ -100,6 +100,6 @@
 !============================================================================
 !============================================================================
-subroutine updaterice_microco2(qice, qccn, nccn, temperature, coeff, rice, rhocloudco2)
-use tracer_mod, only: rho_dust
+subroutine updaterice_microco2(qice, qccn, nccn, qccn_h2o, nccn_h2o, temperature, coeff, rice, rhocloudco2)
+use tracer_mod, only: rho_dust, rho_ice
 use comcstfi_h, only: pi
 use density_co2_ice_mod, only: density_co2_ice
@@ -108,10 +108,10 @@
 !CO2 clouds parameter update by CL and JA 09/16
 
-DOUBLE PRECISION, intent(in)  :: qice,qccn,nccn
+DOUBLE PRECISION, intent(in)  :: qice,qccn,nccn, qccn_h2o, nccn_h2o
 real, intent(in)  :: coeff         ! this coeff is tauscaling if microphy = T (possibly ccn_factor^-1 otherwise)
 real, intent(in)  :: temperature ! temperature effective for density co2_ice computation
 real, intent(out) :: rhocloudco2 ! rhocloud is needed for sedimentation and is also a good diagnostic variable
 double precision, intent(out) :: rice
-real nccn_true,qccn_true ! nombre et masse de CCN 
+real nccn_true,qccn_true, qccn_h2o_true, nccn_h2o_true ! nombre et masse de CCN
 double precision :: rho_ice_co2T ! density co2_ice Temperature-dependent
 
@@ -119,17 +119,20 @@
 qccn_true = max(qccn * coeff, 1.e-30)
 
+nccn_h2o_true = max(nccn_h2o, 1.e-30)
+qccn_h2o_true = max(qccn_h2o, 1.e-30)
+
 call density_co2_ice(dble(temperature), rho_ice_co2T)
 
-rhocloudco2 = (qice * rho_ice_co2T + qccn_true*rho_dust) / (qice + qccn_true)
+rhocloudco2 = (qice * rho_ice_co2T + qccn_true*rho_dust + qccn_h2o_true*rho_ice) / (qice + qccn_true + qccn_h2o_true)
 
 rhocloudco2 = min(max(rhocloudco2,rho_ice_co2T), rho_dust)
 
-rice = (qice + qccn_true) * 0.75 / pi / rhocloudco2 / nccn_true
-
-rhocloudco2 = (qice * rho_ice_co2T + qccn_true*rho_dust) / (qice + qccn_true)
-
-if (rice .le. r3iceco2min) then !r3icemin radius power 3 ?
+rice = (qice + qccn_true + qccn_h2o_true) * 0.75 / pi / rhocloudco2 / (nccn_true + nccn_h2o_true)
+
+rhocloudco2 = (qice * rho_ice_co2T + qccn_true*rho_dust + qccn_h2o_true*rho_ice) / (qice + qccn_true + qccn_h2o_true)
+
+if (rice <= r3iceco2min) then !r3icemin radius power 3 ?
   rice = riceco2min
-else if (rice .ge. r3iceco2max) then !r3icemin radius power 3 ?
+else if (rice >= r3iceco2max) then !r3icemin radius power 3 ?
   rice = riceco2max
 else 
Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 2561)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 2562)
@@ -15,4 +15,7 @@
      &                       igcm_co2_ice, nuiceco2_ref,
      &                       igcm_ccnco2_number, igcm_ccnco2_mass,
+     &                       igcm_ccnco2_h2o_number,
+     &                       igcm_ccnco2_h2o_mass_ice,
+     &                       igcm_ccnco2_h2o_mass_ccn,
      &                       igcm_ccn_number, nuice_ref, varian,
      &                       ref_r0, igcm_dust_submicron,
@@ -96,5 +99,6 @@
 
       LOGICAL,SAVE :: firstcall=.true.
-
+      REAL Nccnco2, Qccnco2
+      REAL Nccnco2_h2o, Qccnco2_h2o
       REAL CBRT
       EXTERNAL CBRT
@@ -191,7 +195,20 @@
           DO l=1,nlayer
             DO ig=1,ngrid
+              Nccnco2 = pq(ig,l,igcm_ccnco2_number)
+              Qccnco2 = pq(ig,l,igcm_ccnco2_mass)
+              Nccnco2_h2o = 0.
+              Qccnco2_h2o = 0.
+              if (co2useh2o) then
+                Nccnco2_h2o = pq(ig,l,igcm_ccnco2_h2o_number)
+                Qccnco2_h2o = pq(ig,l,igcm_ccnco2_h2o_mass_ice)
+     &                        + pq(ig,l,igcm_ccnco2_h2o_mass_ccn)
+                Nccnco2 = Nccnco2 - Nccnco2_h2o
+                Qccnco2 = Qccnco2 - Qccnco2_h2o
+              end if
+
               call updaterice_microco2(dble(pq(ig,l,igcm_co2_ice)),
-     &                              dble(pq(ig,l,igcm_ccnco2_mass)),
-     &                              dble(pq(ig,l,igcm_ccnco2_number)),
+     &                              dble(Qccnco2), dble(Nccnco2),
+     &                              dble(Qccnco2_h2o),
+     &                              dble(Nccnco2_h2o),
      &                              pt(ig,l),
      &                              tauscaling(ig),riceco2(ig,l),
