Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 628)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 629)
@@ -2,5 +2,5 @@
      &                rdust,rice,nuice,
      &                reffrad,nueffrad,
-     &                pq)
+     &                pq,tauscaling)
 
        IMPLICIT NONE
@@ -68,8 +68,13 @@
 
       REAL, PARAMETER :: ccn0 = 1.3E8
-
-c      LOGICAL firstcall
-c      DATA firstcall/.true./
-c      SAVE firstcall
+      
+c     For microphysics only:      
+      REAL Mo,No                       ! Mass and number of ccn
+      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
+      REAL tauscaling(ngridmx)         ! Convertion factor for qccn and Nccn
+
+      LOGICAL firstcall
+      DATA firstcall/.true./
+      SAVE firstcall
 
       REAL CBRT
@@ -81,12 +86,8 @@
 c     ---------------------
 
-c==================================================================
-
-c      IF (firstcall) THEN
-c       At firstcall, rdust and rice are not known; therefore
-c         they need to be computed below.
-
-c      Correction TN 17/04: rdust and rice must be updated at all steps, 
-c      otherwise it is a possible source of bugs
+
+c==================================================================
+c 1. Update radius from fields from dynamics or initial state
+c==================================================================
 
 c       1.1 Dust particles
@@ -110,20 +111,41 @@
           ENDDO
         ENDIF
+        
 c       1.2 Water-ice particles
 c       -----------------------
-        IF (water.AND.activice) THEN
-          DO l=1,nlayer
-            DO ig=1,ngrid
-              rice(ig,l) = max( CBRT(
-     &          (pq(ig,l,igcm_h2o_ice)/rho_ice +
-     &          ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
-     &          (ccn0*4./3.*pi)),rdust(ig,l) )
-              nuice(ig,l) = nuice_ref
-            ENDDO
-          ENDDO
+        IF (water.AND.activice) THEN          
+          IF ((firstcall).or.(microphys.eq..false.)) THEN
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                rice(ig,l) = max(CBRT(
+     &            (pq(ig,l,igcm_h2o_ice)/rho_ice +
+     &            ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
+     &            (ccn0*4./3.*pi)),rdust(ig,l) )
+                nuice(ig,l) = nuice_ref
+              ENDDO
+            ENDDO
+          firstcall = .false.
+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.
+c    Therefore, we use a typical value ccn0 at firstcall, like it is done without microphysics.
+          ELSE
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                Mo = pq(ig,l,igcm_h2o_ice) + 
+     &              pq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
+                No = pq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
+                rhocloud(ig,l) =  pq(ig,l,igcm_h2o_ice)*rho_ice / Mo
+     &           + pq(ig,l,igcm_ccn_mass)*tauscaling(ig)*rho_dust/Mo
+                rhocloud(ig,l) =
+     &            min(max(rhocloud(ig,l),rho_ice),rho_dust)
+                rice(ig,l) =
+     &           CBRT( Mo/No * 0.75 / pi / rhocloud(ig,l))
+                nuice(ig,l) = nuice_ref
+              ENDDO
+            ENDDO
+          ENDIF ! of if ((firstcall).or.(microphys.eq.false))
         ENDIF ! of if (water.AND.activice)
         
-c        firstcall = .false.
-c      ENDIF ! of if firstcall
 
 c==================================================================