Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 532)
+++ trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 541)
@@ -192,4 +192,9 @@
         write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
         print*,'-----------------------------------'
+
+        do i=1,nbin_cld+1
+            rb_cld(i) = dlog(rb_cld(i))  !! we save that so that it is not computed
+                                         !! at each timestep and gridpoint
+        enddo
 
 c       Contact parameter of water ice on dust ( m=cos(theta) )
@@ -265,4 +270,5 @@
 
       Cste = ptimestep * 4. * pi * rho_ice
+      dev2 = 1. / ( sqrt(2.) * sigma_ice )
 
       call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
@@ -286,17 +292,14 @@
         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
         No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
-        dev2 = 1. / ( sqrt(2.) * sigma_ice )
         Rn = rdust(ig,l)
         Rn = -dlog(Rn) 
         Rm = Rn - 3. * sigma_ice*sigma_ice  
-        yeah = dlog(rb_cld(1))
-        yeahn = derf( (yeah+Rn) *dev2)
-        yeahm = derf( (yeah+Rm) *dev2)
+        yeahn = derf( (rb_cld(1)+Rn) *dev2)
+        yeahm = derf( (rb_cld(1)+Rm) *dev2)
         do i = 1, nbin_cld
           n_aer(i) = -0.5 * No * yeahn !! this ith previously computed
           m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed
-          yeah = dlog(rb_cld(i+1))     !! this (i+1)th now computed
-          yeahn = derf( (yeah+Rn) *dev2)
-          yeahm = derf( (yeah+Rm) *dev2)
+          yeahn = derf( (rb_cld(i+1)+Rn) *dev2)
+          yeahm = derf( (rb_cld(i+1)+Rm) *dev2)
           n_aer(i) = n_aer(i) + 0.5 * No * yeahn
           m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm
@@ -363,8 +366,11 @@
      &                      / Mo * rho_dust
           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
-          rice(ig,l) =
-     &      CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
+          if ((Mo.lt.1.e-20) .or. (No.le.1)) then
+              rice(ig,l) = 1.e-8
+          else
+              rice(ig,l) =
+     &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
+          endif
 c          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
-          if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
           seq  = exp( 2.*sig(zt(ig,l))*mh2o / 
      &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
@@ -483,7 +489,10 @@
           
           rice_out(ig,l)=rice(ig,l)
-          rice(ig,l) =
-     &      CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
-          if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
+          if ((Mo.lt.1.e-20) .or. (No.le.1)) then
+              rice(ig,l) = 1.e-8
+          else
+              rice(ig,l) =
+     &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
+          endif
           rice_out(ig,l)=rice(ig,l)-rice_out(ig,l)
           
Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 532)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 541)
@@ -402,6 +402,10 @@
      &                   *tauscaling(ig)/ Mo * rho_dust
           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
-          rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
-           if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8
+          if ((Mo.lt.1.e-15) .or. (No.le.50)) then
+              rice(ig,l) = 1.e-8
+          else
+              !! AS: only perform computations if conditions not met
+              rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
+          endif
          ENDDO
         ENDDO
@@ -414,5 +418,5 @@
             rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
      &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
-     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
+     &      +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l))
      &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
           ENDDO