Index: trunk/LMDZ.MARS/libf/phymars/callradite.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 740)
@@ -389,5 +389,5 @@
      &                rdust,rice,nuice,
      &                reffrad,nueffrad,
-     &                pq,tauscaling)
+     &                pq,tauscaling,tau,pplay)
 
 c     Computing 3D scattering parameters:
Index: trunk/LMDZ.MARS/libf/phymars/callsedim.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 740)
@@ -3,7 +3,8 @@
      &                rsedcloud,rhocloud,
      &                pq, pdqfi, pdqsed,pdqs_sed,nq, 
-     &                tau, tauscaling)
+     &                tau,tauscaling)
 ! to use  'getin'
-      USE ioipsl_getincom 
+      USE ioipsl_getincom
+      USE updaterad
       IMPLICIT NONE
 
@@ -54,7 +55,4 @@
 c   ------
 
-      REAL CBRT
-      EXTERNAL CBRT
-
       INTEGER l,ig, iq
 
@@ -67,6 +65,6 @@
       real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
                                     !   dust (m)
-      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
-                                    !   CCNs (m)
+!      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
+!                                    !   CCNs (m)
 c     Sedimentation radius of water ice
       real rsedcloud(ngridmx,nlayermx)
@@ -186,5 +184,5 @@
         ENDIF !of if (doubleq)
 
-        IF (scavenging) THEN
+        IF (microphys) THEN
           iccn_mass=0
           iccn_number=0
@@ -201,8 +199,8 @@
             write(*,*) 'callsedim: error! could not identify'
             write(*,*) ' tracers for ccn mass and number mixing'
-            write(*,*) ' ratio and scavenging is activated!'
+            write(*,*) ' ratio and microphys is activated!'
             stop
           endif
-        ENDIF !of if (scavenging)
+        ENDIF !of if (microphys)
 
         IF (water) THEN
@@ -255,18 +253,13 @@
         do l=1,nlay
           do ig=1, ngrid
-            r0dust(ig,l) =
-     &        CBRT(r3n_q*zqi(ig,l,idust_mass)/
-     &        max(zqi(ig,l,idust_number),0.01))
-            r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
+     
+         call updaterdust(zqi(ig,l,igcm_dust_mass),
+     &                    zqi(ig,l,igcm_dust_number),r0dust(ig,l),
+     &                    tauscaling(ig))
+          
           end do
         end do
       endif
-      if (scavenging) then
-        do l=1,nlay
-          do ig=1, ngrid
-            r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
-          end do
-        end do
-      endif
+
 
 c =================================================================
@@ -279,9 +272,6 @@
           if ((doubleq.and.
      &        ((iq.eq.idust_mass).or.
-     &         (iq.eq.idust_number)))) then !.or.
-!     &        (scavenging.and.
-!     &        ((iq.eq.iccn_mass).or.
-!     &        (iq.eq.iccn_number)))) then
-
+     &         (iq.eq.idust_number)))) then
+     
 c           Computing size distribution:
 c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -294,12 +284,4 @@
               end do
               sigma0 = varian
-c            else             ! ccn
-c              do  l=1,nlay
-c                do ig=1, ngrid
-c                  r0(ig,l)=r0ccn(ig,l)
-c                end do
-c              end do
-c              sigma0 = sqrt(log(1.+nuice_sed))
-c            endif
 
 c        Computing mass mixing ratio for each particle size
@@ -358,15 +340,8 @@
 
           do ir=1,nr
-!             IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
-              ! Dust sedimentation
+          
                call newsedim(ngrid,nlay,1,1,ptimestep,
      &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
      &         wq,0.5)
-!             ELSE
-!               ! CCN sedimentation
-!               call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
-!     &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
-!     &         wq,beta)
-!             ENDIF
 
 c            Tendencies
@@ -385,5 +360,4 @@
 c         WATER CYCLE CASE
 c -----------------------------------------------------------------
-c          else if (water.and.(iq.eq.igcm_h2o_ice)) then
            else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
      &       .or. (iq .eq. igcm_h2o_ice)) then
@@ -437,8 +411,11 @@
        DO l = 1, nlay
         DO ig=1,ngrid
-          rdust(ig,l)=
-     &      CBRT(r3n_q*zqi(ig,l,idust_mass)/
-     &      max(zqi(ig,l,idust_number),0.01))
-          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
+        
+     
+         call updaterdust(zqi(ig,l,igcm_dust_mass),
+     &                    zqi(ig,l,igcm_dust_number),rdust(ig,l),
+     &                    tauscaling(ig))     
+
+          
         ENDDO
        ENDDO
@@ -448,35 +425,28 @@
 c     -------------------------------------
       if (water) then
-       IF(scavenging) THEN 
+       IF(microphys) THEN 
+       
+       
         DO l = 1, nlay
           DO ig=1,ngrid
-            Mo   = zqi(ig,l,igcm_h2o_ice) + 
-     &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
-            No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
-            rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
-     &                       +zqi(ig,l,iccn_mass)* 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.1)) rice(ig,l) = 1.e-8
-!           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
-!           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
+
+         call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
+     &    zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
+     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
            
           ENDDO
         ENDDO
+        
        ELSE
+       
         DO l = 1, nlay
           DO ig=1,ngrid
-            ccntyp = 
-     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
-            ccntyp = ccntyp /ccn_factor
-            rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
-     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
-     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
+          
+            call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
+     &                      tau(ig,1),zlay(ig,l),rice(ig,l))
+
           ENDDO
         ENDDO
-       ENDIF ! of IF(scavenging)
+       ENDIF ! of IF(microphys)
       endif ! of if (water)
 
Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 740)
@@ -5,4 +5,5 @@
 ! to use  'getin'
       USE ioipsl_getincom
+      USE updaterad
       implicit none
       
@@ -69,11 +70,7 @@
       REAL*8   derf ! Error function
       !external derf
-
-      REAL CBRT
-      EXTERNAL CBRT
-
+     
       INTEGER ig,l,i
       
-
       REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
       REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
@@ -104,5 +101,5 @@
 
       REAL res      ! Resistance growth
-
+      
 
 c     Parameters of the size discretization
@@ -131,5 +128,9 @@
       
       LOGICAL test_flag    ! flag for test/debuging outputs
-      SAVE    test_flag     
+      SAVE    test_flag    
+
+
+      REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) 
+      REAL res_out(ngrid,nlay)
  
 
@@ -185,5 +186,5 @@
 
         do i=1,nbin_cld+1
-    !        rb_cld(i) = log(rb_cld(i))  
+!           rb_cld(i) = log(rb_cld(i))  
             rb_cld(i) = dlog(rb_cld(i))  !! we save that so that it is not computed
                                          !! at each timestep and gridpoint
@@ -217,4 +218,7 @@
       cste = 4*pi*rho_ice*ptimestep
 
+      res_out(:,:) = 0
+      rice(:,:) = 1.e-8
+
 c     Initialize the tendencies
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
@@ -268,9 +272,6 @@
        IF ( satu .ge. 1. ) THEN         ! if there is condensation
 
-c         Update rdust from last tendencies
-        rdust(ig,l)=
-     &     CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
-     &      (zq(ig,l,igcm_dust_number)+1.e-30))
-        rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
+        call updaterccn(zq(ig,l,igcm_dust_mass),
+     &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
 
 
@@ -332,12 +333,12 @@
 c       Update Dust particles
         zq(ig,l,igcm_dust_mass)   = 
-     &         zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig)
+     &  zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 
         zq(ig,l,igcm_dust_number) = 
-     &         zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
+     &  zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
 c       Update CCNs
         zq(ig,l,igcm_ccn_mass)   = 
-     &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
+     &  zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
         zq(ig,l,igcm_ccn_number) = 
-     &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
+     &  zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
         
         ENDIF ! of is satu >1
@@ -353,29 +354,21 @@
        IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth 
 
-
-        Mo   = zq(ig,l,igcm_h2o_ice) + 
-     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
+     
+        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))
+
         No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
 
-
-        rhocloud(ig,l) =  zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
-     &                  + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
-     &                  / Mo * rho_dust
-        rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
-          
-
-c       ice crystal radius    
-        rice (ig,l) = 
-     &      CBRT( real(Mo/No) * 0.75 / pi / rhocloud(ig,l) )
-
-
 c       saturation at equilibrium
-        seq  = exp( 2.*sig(zt(ig,l))*mh2o / 
-     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
-
+c       rice should not be too small, otherwise seq value is not valid
+        seq  = exp(2.*sig(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*
+     &             max(rice(ig,l),1.e-7)))
+        
 c       get resistance growth
         call growthrate(zt(ig,l),pplay(ig,l),
      &             real(ph2o/satu),rice(ig,l),res)
 
+        res_out(ig,l) = res
 
 ccccccc  implicit scheme of mass growth
@@ -434,20 +427,43 @@
       
       ! Get cloud tendencies
-      pdqcloud(1:ngrid,1:nlay,1:nq) =
-     & (zq(1:ngrid,1:nlay,1:nq)-zq0(1:ngrid,1:nlay,1:nq))/ptimestep
-     
-      
+        pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
+     &   (zq(1:ngrid,1:nlay,igcm_h2o_vap) - 
+     &    zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
+        pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
+     &   (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
+     &    zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
+        pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
+     &   (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
+     &    zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
+        pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
+     &   (zq(1:ngrid,1:nlay,igcm_ccn_number) -
+     &    zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
+     
+      if (scavenging) then
+      
+        pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
+     &   (zq(1:ngrid,1:nlay,igcm_dust_mass) -
+     &    zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
+        pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
+     &   (zq(1:ngrid,1:nlay,igcm_dust_number) -
+     &    zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
+          
+      endif
+     
+     
      
 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
-!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
+!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
       IF (test_flag) then
       
 !       error2d(:) = 0.
-!       DO l=1,nlay
-!       DO ig=1,ngrid
+       DO l=1,nlay
+       DO ig=1,ngrid
 !         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
-!       ENDDO
-!       ENDDO
+          satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
+          satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
+       ENDDO
+       ENDDO
 
       print*, 'count is ',countcells, ' i.e. ',
@@ -470,20 +486,22 @@
 !         call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
 !     &                    error_out)
-!         call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
-!     &                    satubf)
-!         call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
-!     &                    satuaf)
-!         call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
-!     &                    zq0(1,:,igcm_h2o_vap))
-!         call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
-!     &                    zq(1,:,igcm_h2o_vap))
-!         call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
-!     &                    zq0(1,:,igcm_h2o_ice))
-!         call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
-!     &                    zq(1,:,igcm_h2o_ice))
-!         call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
-!     &                    ccnbf)
-!         call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
-!     &                    ccnaf)
+         call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
+     &                    res_out)
+         call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
+     &                    satubf)
+         call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
+     &                    satuaf)
+         call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
+     &                    zq0(1,:,igcm_h2o_vap))
+         call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
+     &                    zq(1,:,igcm_h2o_vap))
+         call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
+     &                    zq0(1,:,igcm_h2o_ice))
+         call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
+     &                    zq(1,:,igcm_h2o_ice))
+         call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
+     &                    zq0(1,:,igcm_ccn_number))
+         call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
+     &                    zq(1,:,igcm_ccn_number))
 c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
 c     &                    gr_out)
@@ -494,10 +512,10 @@
 c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
 c     &                    dN_out)
-!         call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
-!     &                    zqsat)
+         call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
+     &                    zqsat)
 !         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
 !     &                    satu_out)
-!         call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
-!     &                    rice)
+         call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
+     &                    rice)
 !         call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
 !     &                    rdust)
Index: trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 740)
@@ -477,23 +477,23 @@
 ! if scavenging is used, then dustbin should be > 0
 
-         if (scavenging.and.(dustbin.lt.1)) then
-           print*,'if scavenging is used, then dustbin should > 0'
-           stop
-         endif
          if ((microphys.and..not.doubleq).or.
-     &       (microphys.and..not.active).or.
-     &       (microphys.and..not.scavenging).or.
      &       (microphys.and..not.water)) then
-           print*,'if microphys is used, then doubleq,'
-           print*,'active, water and scavenging must all be used!'
-           stop
-         endif
-         if ((scavenging.and..not.doubleq).or.
-     &       (scavenging.and..not.active).or.
-     &       (scavenging.and..not.water).or.
-     &       (scavenging.and..not.microphys)) then
-           print*,'if scavenging is used, then doubleq,'
-           print*,'active, water and microphys must all be used!'
-           stop
+             print*,'if microphys is used, then doubleq,'
+             print*,'and water must be used!'
+             stop
+         endif
+         if (microphys.and..not.scavenging) then
+             print*,''
+             print*,'----------------WARNING-----------------'
+             print*,'microphys is used without scavenging !!!'
+             print*,'----------------WARNING-----------------'
+             print*,''
+         endif
+
+         if ((scavenging.and..not.microphys).or.
+     &       (scavenging.and.(dustbin.lt.1))) then
+             print*,'if scavenging is used, then microphys'
+             print*,'must be used!'
+             stop
          endif
 
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 740)
@@ -162,5 +162,5 @@
         enddo
       endif ! of if (doubleq)
-      if (scavenging) then
+      if (microphys) then
         do iq=1,nqmx
           if (noms(iq).eq."ccn_mass") then
@@ -173,5 +173,5 @@
           endif
         enddo
-      endif ! of if (scavenging)
+      endif ! of if (microphys)
       if (submicron) then
         do iq=1,nqmx
Index: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 740)
@@ -3,4 +3,5 @@
      &             pq,pdq,pdqcloud,pdtcloud,
      &             nq,tau,rice)
+      USE updaterad
       implicit none
 c------------------------------------------------------------------
@@ -41,5 +42,4 @@
       integer nq                 ! nombre de traceurs
       REAL ptimestep             ! pas de temps physique (s)
-!      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
       REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
@@ -68,6 +68,6 @@
       SAVE firstcall
 
-      REAL CBRT
-      EXTERNAL CBRT
+            
+      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
 
       INTEGER ig,l
@@ -86,5 +86,4 @@
 c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
       
-      REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
 
 c-----------------------------------------------------------------------
@@ -110,31 +109,4 @@
       enddo
 
-      do l=1,nlay
-        do ig=1,ngrid
-
-c         Typical dust radius profile:
-          rdusttyp(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
-
-c         Typical CCN profile:
-c         Corrected equation, following Montmessin et al. 2004
-c           (N0=2e6 m-3 has been converted into N0=1.3e8 kg-1, otherwise
-c            the equation for rice is not homogeneous...)
-          ccntyp(ig,l)=
-     &       1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
-c         The previously used profile was not correct:
-c         ccntyp(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) *
-c    &       2.e+6/0.1*max(tau(ig,1),0.001)*exp(-pzlay(ig,l)/10000.)
-
-c         Reduce number of nuclei
-!         TEMPORAIRE : decrease the number of CCNs FF 04/200
-!         reduction facteur 3
-!         ccntyp(ig,l) = ccntyp(ig,l) / 27.
-!         reduction facteur 2
-!         ccntyp(ig,l) = ccntyp(ig,l) / 8.
-c -----------------------------------------------------------------
-          ccntyp(ig,l) = ccntyp(ig,l) / ccn_factor
-
-        enddo ! of do ig=1,ngrid
-      enddo ! of do l=1,nlay
 
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
@@ -167,15 +139,4 @@
           zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
           
-          Mcon_out(ig,l) = dzq
-
-c         Calcul du rayon moyen des particules de glace.
-c         Hypothese : Dans une couche, la glace presente se
-c         repartit uniformement autour du nbre de poussieres
-c         dont le rayon moyen est prescrit par rdusttyp.
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          rice(ig,l)=max( CBRT ( (zq(ig,l,igcm_h2o_ice)/rho_ice
-     &      +ccntyp(ig,l)*(4./3.)*pi*rdusttyp(ig,l)**3.)
-     &      /(ccntyp(ig,l)*4./3.*pi) ), rdusttyp(ig,l)) 
-
 
         enddo ! of do ig=1,ngrid
@@ -195,14 +156,12 @@
       end do
 
-c------------------------------------------------------------------
-c     TEST_JBM
-!      IF (ngrid.eq.1) THEN
-!         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
-!     &                    Mcon_out)
-!         call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1,
-!     &                    rdusttyp)
-!         call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1,
-!     &                    ccntyp)
-!      ENDIF
+c     ice crystal radius
+      do l=1, nlay
+        do ig=1,ngridmx
+          call updaterice_typ(zq(ig,l,igcm_h2o_ice),
+     &       tau(ig,1),pzlay(ig,l),rice(ig,l))
+        end do
+      end do
+
 c------------------------------------------------------------------
       return
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 740)
@@ -307,12 +307,10 @@
           longwatercaptag(203)    = .true. ! outlier, lat = 75
           longwatercaptag(207)    = .true. ! outlier, lat = 75
-          longwatercaptag(242)    = .true. ! outlier, lat = 75
           longwatercaptag(244)    = .true. ! outlier, lat = 75
           longwatercaptag(246)    = .true. ! outlier, lat = 75
-          longwatercaptag(248)    = .true. ! outlier, lat = 75
           longwatercaptag(250)    = .true. ! outlier, lat = 75
           longwatercaptag(252)    = .true. ! outlier, lat = 75
           longwatercaptag(254)    = .true. ! outlier, lat = 75
-          longwatercaptag(256:257)= .true. ! outlier, lat = 75
+          longwatercaptag(256)    = .true. ! outlier, lat = 75
           endif
 !--------------------------------------------------------------       
Index: trunk/LMDZ.MARS/libf/phymars/updaterad.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updaterad.F90	(revision 740)
+++ trunk/LMDZ.MARS/libf/phymars/updaterad.F90	(revision 740)
@@ -0,0 +1,279 @@
+module updaterad
+
+
+
+! This module intents to group together all ice and dust radius computations done in the GCM,
+! so that it stays coherent through the code and make (numerous) radius bugs easier to debug.
+! All different thresholds values are defined below.
+
+! Radius computation do not always occur on the whole grid (cf. improvedcloud).
+! So, subroutines are designed for scalar values instead of tables
+
+
+! T. Navarro, June 2012
+
+
+
+! For instance, if R^3 is lower than r3icemin, then R is set to ricemin.
+! So, ricemin should be the cubic root of r3icemin, but not necessarily ...
+
+real, parameter :: r3icemin = 1.e-30   ! ie ricemin  = 0.0001 microns
+real, parameter :: ricemin  = 1.e-10
+
+real, parameter :: r3icemax = 125.e-12 ! ie ricemax  = 500 microns
+real, parameter :: ricemax  = 500.e-6 
+
+real, parameter :: qice_threshold  = 1.e-15 ! 1.e-10
+
+
+
+real, parameter :: nccn_threshold  =  1.
+real, parameter :: qccn_threshold  =  1.e-20
+
+real, parameter :: r3ccnmin = 1.e-21    ! ie rccnmin = 0.1 microns
+real, parameter :: rccnmin  = 0.1e-6
+
+real, parameter :: r3ccnmax = 125.e-12  ! ie rccnmax  = 500 microns
+real, parameter :: rccnmax  = 500.e-6  
+
+
+
+real, parameter :: ndust_threshold  = 1.
+real, parameter :: qdust_threshold  = 1.e-20
+
+real, parameter :: r3dustmin = 1.e-21  ! ie rdustmin = 0.1 microns
+real, parameter :: rdustmin  = 0.1e-6
+
+real, parameter :: r3dustmax = 125.e-12 ! ie rdustmax  = 500 microns
+real, parameter :: rdustmax  = 500.e-6  
+
+
+real, parameter :: rdust0 = 0.8e-7
+
+
+      
+
+contains
+
+
+!============================================================================
+!============================================================================
+!============================================================================
+! Update ice radius if microphys == true 
+subroutine updaterice_micro(qice,qccn,nccn,coeff,rice,rhocloud)
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "tracer.h"
+
+real, intent(in)  :: qice,qccn,nccn
+real, intent(in)  :: coeff         ! this coeff is tauscaling if microphy = T (possibly ccn_factor^-1 otherwise)
+real, intent(out) :: rice,rhocloud ! rhocloud is needed for sedimentation and is also a good diagnostic variable
+real nccn_true,qccn_true ! radius of the ice crystal core to the power of 3
+
+    
+nccn_true = max(nccn * coeff, 1.e-30)
+qccn_true = max(qccn * coeff, 1.e-30)
+
+
+!! Nota: It is very dangerous to apply a threshold on qccn or nccn to force rice to be ricemin.
+!! Indeed, one can obtain ricemin for small but non negligible qice values, and therefore hugely opaque clouds.
+ 
+if (qice .ge. qice_threshold) then
+
+  rhocloud = (qice*rho_ice + qccn_true*rho_dust) / (qice + qccn_true)
+  rhocloud = min(max(rhocloud,rho_ice),rho_dust)
+
+  rice = (qice + qccn_true) * 0.75 / pi / rhocloud / nccn_true
+  
+  if (rice .le. r3icemin) then
+    rice = ricemin
+  else if (rice .ge. r3icemax) then
+    rice = ricemax
+  else
+    rice = rice**(1./3.) ! here rice is always positive
+  endif
+
+else
+
+  rice = ricemin
+  rhocloud = rho_dust
+
+endif
+
+end subroutine updaterice_micro
+!============================================================================
+!============================================================================
+!============================================================================
+
+
+
+
+
+!============================================================================
+!============================================================================
+!============================================================================
+! Update ice radius from a typical profile if microphys == false
+subroutine updaterice_typ(qice,tau,pzlay,rice)
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "tracer.h"
+
+real, intent(in)  :: qice
+real, intent(in)  :: tau   ! tau for dust
+real, intent(in)  :: pzlay ! altitude at the middle of the layers
+real, intent(out) :: rice
+real rccn,nccn ! radius  and number of ice crystals
+
+!   Typical CCN profile following Montmessin et al. 2004
+!   (N0=2e6 m-3 has been converted into N0=1.3e8 kg-1, otherwise the equation for rice is not homogeneous...)
+     nccn  = 1.3e+8*max(tau,0.001)/0.1*exp(-pzlay/10000.)
+!   The previously used profile was not correct:
+!   Nccn=( epaisseur/masse ) * 2.e+6/0.1*max(tau,0.001)*exp(-pzlay/10000.)
+     
+if (nccn .le. 1) then
+
+  rice = ricemin
+
+else
+
+! Typical dust radius profile:
+  rccn = max(rdust0*exp(-pzlay/18000.),1.e-9)
+  rice  = qice * 0.75 / pi / rho_ice / nccn + rccn*rccn*rccn
+       
+  if (rice .le. r3icemin) then
+    rice = ricemin
+  else if (rice .ge. r3icemax) then
+    rice = ricemax
+  else
+    rice = rice**(1./3.) ! here rice is always positive
+  endif
+
+endif
+ 
+end subroutine updaterice_typ
+!============================================================================
+!============================================================================
+!============================================================================
+
+
+
+
+
+
+!============================================================================
+!============================================================================
+!============================================================================
+! This subroutine computes the geometric mean radius(or number median radius)
+! For a lognormal distribution :
+! geometric mean radius = mass mean radius x exp(-1.5 sigma0^2)
+! To be used with doubleq == true. otherwise, rdust is constant !!!
+subroutine updaterdust(qdust,ndust,rdust,tauscaling)
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "tracer.h"
+
+real, intent(in) :: qdust,ndust ! needed if doubleq
+real, intent(in), optional :: tauscaling ! useful for realistic thresholds
+real, intent(out) :: rdust
+real coeff
+
+if (present(tauscaling)) then
+  coeff = tauscaling ! thresholds on realistic values
+else
+  coeff = 1. ! thresholds on virtual values
+endif
+
+if ((ndust .le. ndust_threshold/coeff) .or. (qdust .le. qdust_threshold/coeff)) then
+
+   rdust = rdustmin
+
+else
+
+   rdust  = r3n_q * qdust / ndust
+
+   if (rdust .le. r3dustmin) then
+       rdust = rdustmin
+   else if (rdust .ge. r3dustmax) then
+       rdust = rdustmax
+   else
+      rdust = rdust**(1./3.)
+   endif
+
+endif
+     
+
+end subroutine updaterdust
+!============================================================================
+!============================================================================
+!============================================================================
+
+
+
+
+
+
+!============================================================================
+!============================================================================
+!============================================================================
+! This subroutine computes the mass mean radius, 
+! used for heterogenous nucleation on CCNs in microphysics.
+! For a lognormal distribution :
+! geometric mean radius = mass mean radius x exp(-1.5 sigma0^2)
+subroutine updaterccn(qccn,nccn,rccn,tauscaling)
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "tracer.h"
+
+real, intent(in) :: qccn,nccn ! needed if doubleq
+real, intent(in), optional :: tauscaling ! useful for realistic thresholds
+
+real, intent(out) :: rccn
+
+real coeff
+
+
+if (present(tauscaling)) then
+  coeff = tauscaling ! threshold on realistic values
+else
+  coeff = 1. ! threshold on virtual values
+endif
+
+if ((nccn .le. nccn_threshold/coeff) .or. (qccn .le. qccn_threshold/coeff)) then
+
+   rccn = rccnmin
+
+else
+
+   rccn  =  qccn * 0.75 / pi / rho_dust / nccn
+
+   if (rccn .le. r3ccnmin) then
+       rccn = rccnmin
+   else if (rccn .ge. r3ccnmax) then
+       rccn = rccnmax
+   else
+       rccn = rccn**(1./3.)
+   endif
+
+endif
+     
+end subroutine updaterccn
+!============================================================================
+!============================================================================
+!============================================================================
+
+
+
+end module updaterad
+
Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 737)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 740)
@@ -6,4 +6,5 @@
 ! to use  'getin'
       USE ioipsl_getincom
+      USE updaterad
       IMPLICIT NONE
 
@@ -92,7 +93,4 @@
       REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
       REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
-      
-      REAL CBRT
-      EXTERNAL CBRT
 
 
@@ -261,4 +259,5 @@
           ENDDO
        ENDDO
+       
 c------ Tracers tendencies pdqcloud
        DO l=1,nlay
@@ -270,11 +269,10 @@
      &        subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
      &       - pdq(ig,l,igcm_h2o_vap)
-          IF(microphys) THEN
-            pdqcloud(ig,l,igcm_dust_mass) = 
-     &        subpdq(ig,l,igcm_dust_mass)/real(imicro) 
-     &       - pdq(ig,l,igcm_dust_mass)
-            pdqcloud(ig,l,igcm_dust_number) = 
-     &        subpdq(ig,l,igcm_dust_number)/real(imicro) 
-     &       - pdq(ig,l,igcm_dust_number)
+         ENDDO
+       ENDDO
+       
+       IF(microphys) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
             pdqcloud(ig,l,igcm_ccn_mass) = 
      &        subpdq(ig,l,igcm_ccn_mass)/real(imicro) 
@@ -283,12 +281,50 @@
      &        subpdq(ig,l,igcm_ccn_number)/real(imicro) 
      &       - pdq(ig,l,igcm_ccn_number)
-           ENDIF
-         ENDDO
-       ENDDO
+         ENDDO
+        ENDDO
+       ENDIF
+       
+       IF(scavenging) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,igcm_dust_mass) = 
+     &        subpdq(ig,l,igcm_dust_mass)/real(imicro) 
+     &       - pdq(ig,l,igcm_dust_mass)
+            pdqcloud(ig,l,igcm_dust_number) = 
+     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
+     &       - pdq(ig,l,igcm_dust_number)
+         ENDDO
+        ENDDO
+       ENDIF
 
 c------- Due to stepped entry, other processes tendencies can add up to negative values
 c------- Therefore, enforce positive values and conserve mass
 
+
        IF(microphys) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+          IF ((pq(ig,l,igcm_ccn_number) + 
+     &      ptimestep* (pdq(ig,l,igcm_ccn_number) + 
+     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
+     &   .or. (pq(ig,l,igcm_ccn_mass) + 
+     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) + 
+     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
+         pdqcloud(ig,l,igcm_ccn_number) =
+     &     - pq(ig,l,igcm_ccn_number)/ptimestep 
+     &     - pdq(ig,l,igcm_ccn_number) + 1.
+         pdqcloud(ig,l,igcm_dust_number) =  
+     &     -pdqcloud(ig,l,igcm_ccn_number)
+         pdqcloud(ig,l,igcm_ccn_mass) =
+     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
+     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
+         pdqcloud(ig,l,igcm_dust_mass) = 
+     &     -pdqcloud(ig,l,igcm_ccn_mass)
+          ENDIF
+         ENDDO
+        ENDDO
+       ENDIF
+
+       IF(scavenging) THEN
         DO l=1,nlay
          DO ig=1,ngrid
@@ -310,21 +346,4 @@
      &     -pdqcloud(ig,l,igcm_dust_mass)
           ENDIF
-          IF ((pq(ig,l,igcm_ccn_number) + 
-     &      ptimestep* (pdq(ig,l,igcm_ccn_number) + 
-     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
-     &   .or. (pq(ig,l,igcm_ccn_mass) + 
-     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) + 
-     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
-         pdqcloud(ig,l,igcm_ccn_number) =
-     &     - pq(ig,l,igcm_ccn_number)/ptimestep 
-     &     - pdq(ig,l,igcm_ccn_number) + 1.
-         pdqcloud(ig,l,igcm_dust_number) =  
-     &     -pdqcloud(ig,l,igcm_ccn_number)
-         pdqcloud(ig,l,igcm_ccn_mass) =
-     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
-     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
-         pdqcloud(ig,l,igcm_dust_mass) = 
-     &     -pdqcloud(ig,l,igcm_ccn_mass)
-          ENDIF
          ENDDO
         ENDDO
@@ -358,86 +377,56 @@
          DO ig=1,ngrid
 
-          rdust(ig,l)=
-     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
-     &           (pdq(ig,l,igcm_dust_mass)
-     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
-     &      (pq(ig,l,igcm_dust_number)+
-     &           (pdq(ig,l,igcm_dust_number)+
-     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
-          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
-
-          Mo = pq(ig,l,igcm_h2o_ice) 
-     &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
-     &         + (pq(ig,l,igcm_ccn_mass)
-     &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
-     &         *tauscaling(ig) + 1.e-30
-          No = (pq(ig,l,igcm_ccn_number)
-     &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
-     &         *tauscaling(ig) + 1.e-30
-          rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) +  
-     &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
-     &                   / Mo * rho_ice
-     &                  + (pq(ig,l,igcm_ccn_mass)
-     &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
-     &                   *tauscaling(ig)/ Mo * rho_dust
-          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
-          if ((Mo.lt.1.e-15) .or. (No.le.1.)) 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.)
-              rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
-          endif
-          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
-     &                         (1.+nuice_sed)*(1.+nuice_sed),
-     &                         rdust(ig,l) )
-          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
-         ENDDO
-        ENDDO
-       ELSE
-        
-        if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0)) then
-         ! recompute rdust(), if we have dust_mass & dust_number tracers
-         DO l=1,nlay
+        call updaterdust(
+     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
+     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
+     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
+     &    pq(ig,l,igcm_dust_number) +                 ! dust number
+     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
+     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
+     &    rdust(ig,l))
+
+         ENDDO
+        ENDDO
+      ENDIF
+        
+        
+      IF(microphys) THEN
+       
+       DO l=1, nlay
+         DO ig=1,ngrid
+
+        call updaterice_micro(
+     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
+     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
+     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
+     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
+     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
+     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
+     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
+     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+          
+         ENDDO
+       ENDDO
+        
+      ELSE ! no microphys
+        
+        DO l=1,nlay
           DO ig=1,ngrid
-           rdust(ig,l)=
-     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
-     &           (pdq(ig,l,igcm_dust_mass)
-     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
-     &      (pq(ig,l,igcm_dust_number)+
-     &           (pdq(ig,l,igcm_dust_number)+
-     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
-           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
-          ENDDO
-         ENDDO
-        endif ! of if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0))
-        
-        DO l=1,nlay
-          DO ig=1,ngrid
-            ccntyp = 
-     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
-            ccntyp = ccntyp /ccn_factor
-            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)*rdust(ig,l)*rdust(ig,l))
-     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
-            rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
-     &                         (1.+nuice_sed)*(1.+nuice_sed),
-     &                         rdust(ig,l) )
-            rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
-          ENDDO
-        ENDDO
-       ENDIF ! of IF(scavenging)
-       
-       
-! used for rad. transfer calculations
-! nuice is constant because a lognormal distribution is prescribed
-      nuice(1:ngrid,1:nlay)=nuice_ref 
-
-
-!--------------------------------------------------------------
-!--------------------------------------------------------------
-      
-
+          
+        call updaterice_typ(
+     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
+     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &    tau(ig,l),pzlay(ig,l),rice(ig,l)) 
+
+          ENDDO
+         ENDDO
+       
+       ENDIF ! of IF(microphys)
+      
+      
+      
 c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -450,4 +439,20 @@
         end if
       end do
+       
+       
+       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))
+!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
+         ENDDO
+       ENDDO
+       
+! used for rad. transfer calculations
+! nuice is constant because a lognormal distribution is prescribed
+      nuice(1:ngrid,1:nlay)=nuice_ref 
+
+
 
 c=======================================================================
