Index: trunk/LMDZ.MARS/libf/phymars/callsedim.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 626)
@@ -206,15 +206,15 @@
         ENDIF !of if (scavenging)
 
-!        IF (water) THEN
-!         write(*,*) "correction for the shape of the ice particles ?"
-!         beta=0.75 ! default value
-!         call getin("ice_shape",beta)
-!         write(*,*) " ice_shape = ",beta
-!
-!          write(*,*) "water_param nueff Sedimentation:", nuice_sed
-!          IF (activice) THEN
-!            write(*,*) "water_param nueff Radiative:", nuice_ref
-!          ENDIF
-!        ENDIF
+        IF (water) THEN
+         write(*,*) "correction for the shape of the ice particles ?"
+         beta=0.75 ! default value
+         call getin("ice_shape",beta)
+         write(*,*) " ice_shape = ",beta
+
+          write(*,*) "water_param nueff Sedimentation:", nuice_sed
+          IF (activice) THEN
+            write(*,*) "water_param nueff Radiative:", nuice_ref
+          ENDIF
+        ENDIF
       
         firstcall=.false.
@@ -383,29 +383,29 @@
           enddo ! of do ir=1,nr
 c -----------------------------------------------------------------
-c         WATER CYCLE CASE --- NOW DONE IN WATERCLOUD !!
-c -----------------------------------------------------------------
-!          else if (water.and.(iq.eq.igcm_h2o_ice)) then
-!            if (microphys) then
-!              ! water ice sedimentation
-!              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
-!     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
-!     &        zqi(1,1,iq),wq,beta)
-!            else
-!              ! water ice sedimentation
-!              call newsedim(ngrid,nlay,ngrid*nlay,1,
-!     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
-!     &        zqi(1,1,iq),wq,beta)
-!            endif ! of if (microphys)
+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
+            if (microphys) then
+              ! water ice sedimentation
+              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
+     &        zqi(1,1,iq),wq,beta)
+            else
+              ! water ice sedimentation
+              call newsedim(ngrid,nlay,ngrid*nlay,1,
+     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
+     &        zqi(1,1,iq),wq,beta)
+            endif ! of if (microphys)
 c           Tendencies
 c           ~~~~~~~~~~
-!            do ig=1,ngrid 
-!              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
-!            end do
+            do ig=1,ngrid 
+              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
+            end do
 c -----------------------------------------------------------------
 c         GENERAL CASE
 c -----------------------------------------------------------------
-c          else
-           else if ((iq .ne. iccn_mass) .and. (iq .ne. iccn_number)
-     &       .and. (iq .ne. igcm_h2o_ice)) then ! because it is done in watercloud
+          else
             call newsedim(ngrid,nlay,1,1,ptimestep,
      &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
@@ -445,35 +445,35 @@
 c     Update the ice particle size "rice"
 c     -------------------------------------
-!      IF(scavenging) 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.50)) rice(ig,l) = 1.e-8
-c           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
-c           print*, "rice, rho apres", 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))
-!          ENDDO
-!        ENDDO
-!      ENDIF
+      IF(scavenging) 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)
+           
+          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))
+          ENDDO
+        ENDDO
+      ENDIF
 
       RETURN
Index: trunk/LMDZ.MARS/libf/phymars/growthrate.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 626)
@@ -1,4 +1,3 @@
-      subroutine growthrate(timestep,temp,pmid,ph2o,psat,
-     &                      seq,rcrystal,growth)
+      subroutine growthrate(temp,pmid,psat,seq,rcrystal,coeff1,coeff2)
 
       IMPLICIT NONE
@@ -10,4 +9,6 @@
 c     Authors: F. Montmessin
 c       Adapted for the LMD/GCM by J.-B. Madeleine (October 2011)
+c       Use of resistances in the analytical function 
+c            instead of growth rate - T. Navarro (2012)
 c     
 c=======================================================================
@@ -27,12 +28,14 @@
 c   ----------
 
-      REAL timestep
+c     Input
       REAL temp     ! temperature in the middle of the layer (K)
       REAL pmid     ! pressure in the middle of the layer (K)
-      REAL*8 ph2o   ! water vapor partial pressure (Pa)
       REAL*8 psat   ! water vapor saturation pressure (Pa) 
+      REAL*8 seq    ! Equilibrium saturation ratio
       REAL rcrystal ! crystal radius before condensation (m)
-      REAL*8 seq    ! Equilibrium saturation ratio
-      REAL*8 growth 
+
+c     Output
+      REAL coeff1,coeff2     ! coefficients for the analytical relation between time and radius
+
 
 c   local:
@@ -43,4 +46,6 @@
       REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
       REAL*8 Rk,Rd
+      
+      
 
 c-----------------------------------------------------------------------
@@ -73,8 +78,10 @@
      &   ( pi * pmid * (molco2+molh2o)*(molco2+molh2o) 
      &        * sqrt(1.+mh2o/mco2) )
-
+      
       knudsen = temp / pmid * afactor / rcrystal
       lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
-      Dv      = Dv / (1. + lambda * knudsen)
+      
+c      Dv is not corrected. Instead, we use below coefficients coeff1, coeff2
+c      Dv      = Dv / (1. + lambda * knudsen)
 
 c     - Compute Rk
@@ -82,8 +89,13 @@
 c     - Compute Rd
       Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
-
+      
+      
+      coeff1 = real(Rk + Rd*(1. + lambda * knudsen))
+      coeff2 = real(Rk + Rd*(1. - lambda * knudsen))
+      
+c Below are growth rate used for other schemes
 c     - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt)
-      growth = 1. / (seq*Rk+Rd)
-c       growth = (ph2o/psat-seq) / (seq*Rk+Rd)
+c      growth = 1. / (seq*Rk+Rd)
+c      growth = (ph2o/psat-seq) / (seq*Rk+Rd)
 c      rf   = sqrt( max( r**2.+2.*growth*timestep , 0. ) )
 c      dr   = rf-r
Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 626)
@@ -1,4 +1,4 @@
       subroutine improvedclouds(ngrid,nlay,ptimestep,
-     &             pplev,pplay,zlev,pt,pdt,
+     &             pplev,pplay,pt,pdt,
      &             pq,pdq,pdqcloud,pdtcloud,
      &             nq,tauscaling,rdust,rice,nuice,
@@ -21,8 +21,10 @@
 c    values which are then used by the sedimentation and advection
 c    schemes.
+c  A word about the radius growth ...
+c  A word about nucleation and ice growth strategies ...
 
 c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
 c           (October 2011)
-c           T. Navarro, debug & correction (October-December 2011)
+c           T. Navarro, debug,correction, new scheme (October-April 2011)
 c           A. Spiga, optimization (February 2012)
 c------------------------------------------------------------------
@@ -43,5 +45,4 @@
       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
             
       REAL pt(ngrid,nlay)        ! temperature at the middle of the
@@ -95,9 +96,10 @@
       REAL*8 Mo,No
       REAL*8 dN,dM,newvap
-      REAL*8 Rn, Rm, dev2, yeah, yeahn, yeahm
+      REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
       REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
       REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
       REAL*8 rate(nbin_cld)  ! nucleation rate
-      REAL*8 up,dwn,Ctot,gr,seq
+      !REAL*8 up,dwn,Ctot,gr
+      REAl*8 seq
       REAL*8 sig      ! Water-ice/air surface tension  (N.m)
       EXTERNAL sig
@@ -119,28 +121,29 @@
       SAVE sigma_ice
       
+      REAL tdicho, tmax, rmin, rmax, req, rdicho
+      REAL coeff0, coeff1, coeff2
+      REAL error_out(ngridmx,nlayermx)
+      REAL error2d(ngridmx)
+      
+      REAL var1,var2,var3 
+      
+      REAL rccn, epsilon
+      
 c----------------------------------      
-c some outputs for 1D -- TEMPORARY
+c some outputs for 1D -- TESTS
       REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
       REAL dN_out(ngridmx,nlayermx)  ! mass variation for output
       REAL dM_out(ngridmx,nlayermx)  ! number variation for output
-      REAL dM_col(ngridmx)           ! total mass condensed in column
-      REAL dN_col(ngridmx)           ! total mass condensed in column
       REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
-      REAL gr_out(ngridmx,nlayermx)   ! for 1d output
-      REAL newvap_out(ngridmx,nlayermx)   ! for 1d output
-      REAL Mdust_col(ngridmx)        ! total column dust mass
-      REAL Ndust_col(ngridmx)        ! total column dust mass
-      REAL Mccn_col(ngridmx)         ! total column ccn mass
-      REAL Nccn_col(ngridmx)         ! total column ccn mass
-      REAL dMice_col(ngridmx)        ! total column ice mass change
-      REAL drice_col(ngridmx)        ! total column ice radius change
-      REAL icetot(ngridmx)        ! total column ice mass
+      REAL gr_out(ngridmx,nlayermx)   ! for 1d output     
       REAL rice_out(ngridmx,nlayermx) ! ice radius change
       REAL rate_out(ngridmx,nlayermx) ! nucleation rate
+      REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx)
+      REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx)
       INTEGER count
       
-      LOGICAL output_sca     ! scavenging outputs flag for tests
-      
-      output_sca = .false.
+      LOGICAL test_flag    ! flag for test/debuging outputs
+      
+      test_flag = .false.
 c----------------------------------      
 c----------------------------------      
@@ -149,7 +152,7 @@
 
       IF (firstcall) THEN
-
-c       Definition of the size grid
-c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!=============================================================
+! 0. Definition of the size grid
+!=============================================================
 c       rad_cld is the primary radius grid used for microphysics computation.
 c       The grid spacing is computed assuming a constant volume ratio
@@ -215,7 +218,8 @@
       END IF
 
-c-----------------------------------------------------------------------
-c     1. Initialization
-c-----------------------------------------------------------------------
+!=============================================================
+! 1. Initialisation
+!=============================================================
+
 c     Initialize the tendencies
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
@@ -260,15 +264,16 @@
           zq(ig,l,igcm_h2o_ice)=
      &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
-          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 
+          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004 
           zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
         enddo
       enddo
 
-c------------------------------------------------------------------
-c     Cloud microphysical scheme
-c------------------------------------------------------------------
-
-      Cste = ptimestep * 4. * pi * rho_ice
+!=============================================================
+! 2. Compute saturation
+!=============================================================
+
       dev2 = 1. / ( sqrt(2.) * sigma_ice )
+      error_out(:,:) = 0.
+
 
       call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
@@ -283,25 +288,28 @@
         ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
         satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
-        
-        IF ((satu .ge. 1)!        ) THEN             ! if we have condensation
-     &      .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)
-     &             .ge. 2)    ) THEN    ! or sublimation
-
+                
+        IF (( satu .ge. 1. )       .or.                             ! if there is condensation
+     &      ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
+
+
+!=============================================================
+! 3. Nucleation
+!=============================================================
 
 c       Expand the dust moments into a binned distribution
-        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
-        No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
+        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   !+ 1.e-30
+        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
         Rn = rdust(ig,l)
         Rn = -dlog(Rn) 
         Rm = Rn - 3. * sigma_ice*sigma_ice  
-        yeahn = derf( (rb_cld(1)+Rn) *dev2)
-        yeahm = derf( (rb_cld(1)+Rm) *dev2)
+        n_derf = derf( (rb_cld(1)+Rn) *dev2)
+        m_derf = 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
-          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
+          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
+          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
+          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
+          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
+          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
+          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
         enddo
 !!! MORE EFFICIENT COMPUTATIONALLY THAN
@@ -350,5 +358,4 @@
         rate_out(ig,l) = 0
         do i = 1, nbin_cld
-        ! schema simple
           n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
           m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
@@ -358,85 +365,244 @@
         enddo
 
-          Mo   = zq0(ig,l,igcm_h2o_ice) + 
-     &           zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
-          No   = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
-          !write(*,*) "l,cloud particles,cloud mass",l, No, Mo
-          rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
-     &                     +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
-     &                      / Mo * rho_dust
-          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
-          if ((Mo.lt.1.e-20) .or. (No.le.1)) then
-              rice(ig,l) = 1.e-8
+c       Update CCNs, can also be done after the radius growth
+c       Dust particles
+        zq(ig,l,igcm_dust_mass)   = 
+     &         zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig)
+        zq(ig,l,igcm_dust_number) = 
+     &         zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
+c       CCNs
+        zq(ig,l,igcm_ccn_mass)   = 
+     &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
+        zq(ig,l,igcm_ccn_number) = 
+     &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
+        
+
+!=============================================================
+! 4. Ice growth: scheme for radius evolution
+!=============================================================
+
+        Mo   = zq(ig,l,igcm_h2o_ice) + 
+     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig)  + 1.e-30
+        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-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       nuclei radius  
+        rccn  = CBRT(zq(ig,l,igcm_ccn_mass)/
+     &                (pi*rho_dust*zq(ig,l,igcm_ccn_number)*4./3.))
+
+c       ice crystal radius    
+        rice (ig,l) = 
+     &      CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) )
+
+c       enforce physical value : crystal cannot be smaller than its nuclei !
+        rice(ig,l) = max(rice(ig,l), rccn)   
+
+c       saturation at equilibrium
+        seq  = exp( 2.*sig(zt(ig,l))*mh2o / 
+     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
+
+
+c       If there is more than on nuclei, we peform ice growth  
+        var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
+        IF (var1 .ge. -1) THEN
+
+
+      if (test_flag) then
+       print*, ' '
+       print*, ptimestep
+       print*, 'satu,seq', real(satu), real(seq), ig,l
+       print*, 'dN,dM', real(dN),real(dM)           
+       print*,'rccn', rccn
+       print*, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)*tauscaling(ig),
+     &  zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
+      endif
+     
+c         crystal radius to reach saturation at equilibrium (i.e. satu=seq)         
+          req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
+     &          + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)*
+     &           rho_ice/rho_dust - seq * zqsat(ig,l))
+     &         / ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)*
+     &             pi*rho_ice*4./3. )
+          req = CBRT(req)
+          
+c         compute ice radius growth resistances (diffusive and latent heat resistancea)       
+          call growthrate(zt(ig,l),pplay(ig,l),
+     &                    ph2o/satu,seq,req,coeff1,coeff2)
+               
+          coeff0 = -zqsat(ig,l) / (4.*pi*req*rho_ice
+     &               * zq(ig,l,igcm_ccn_number)*tauscaling(ig))
+               
+c         compute tmax, the time needed to reach req            
+          call phi(rice(ig,l),req,coeff1,coeff2,tmax)
+          
+      if (test_flag) then
+          print*, 'coeffs', coeff0,coeff1,coeff2
+          print*, 'req,tmax', req,tmax*coeff0
+          print*, 'i,rmin,rdicho,rmax,tdicho'
+      endif
+
+c         rmin is rice if r increases (satu >1) or req if it decreases (satu<1)
+c         if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn
+          if (satu .ge. seq) then 
+            ! crystal size is increasing
+            rmin = max(min(rice(ig,l),req),rccn) 
+            rmax = max(rice(ig,l),req)
           else
-              rice(ig,l) =
-     &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
+            !  crystal size is decreasing
+            rmax = max(min(rice(ig,l),req),rccn) 
+            rmin = max(rice(ig,l),req)
           endif
-c          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
-          seq  = exp( 2.*sig(zt(ig,l))*mh2o / 
-     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
-
-          call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
-     &                    ph2o,ph2o/satu,seq,rice(ig,l),gr)
-
-          up  = Cste * gr * rice(ig,l) * No * seq + 
-     &            zq(ig,l,igcm_h2o_vap)
-          dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
-  
-          Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
-          
-          newvap = min(up/dwn,Ctot)
-  
-          gr = gr * ( newvap/zqsat(ig,l) - seq )
-
-          
-          dMice = min( max(Cste * No * rice(ig,l) * gr,
-     &                     -zq(ig,l,igcm_h2o_ice) ), 
-     &                 zq(ig,l,igcm_h2o_vap) )
-     
-c----------- TESTS 1D output ---------
-             satu_out(ig,l) = satu
-             Mcon_out(ig,l) = dMice
-             newvap_out(ig,l) = newvap
-             gr_out(ig,l) = gr
-             dN_out(ig,l) = dN
-             dM_out(ig,l) = dM
-c-------------------------------------
-          
-c       Water ice
-          zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) + 
-     &        dMice
-     
-c       Water vapor
-          zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - 
-     &        dMice
-     
+                         !rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm  pour la dicho ...
+          rdicho = 0.5*(rmin+rmax)
+          
+          ! for output
+          var1 = rice(ig,l)
+          var2 = rmin
+          var3 = rmax
+          
+c         Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1        
+          do i = 1,10 ! dichotomy loop
+             
+c            compute tdicho, the time needed to reach rdicho
+             call phi(rdicho,req,coeff1,coeff2,tdicho)
+             !print*, tdicho,tmax
+             tdicho = coeff0*(tdicho - tmax)
+             
+             if (test_flag) print*, i,rmin,rdicho,rmax,tdicho
+                       
+             if (tdicho .ge. ptimestep) then
+                rmax = rdicho
+             else
+                rmin = rdicho
+             endif
+             
+             rdicho = 0.5*(rmin+rmax)
+          
+          enddo  ! of dichotomy loop
+          
+       if (test_flag) then
+        if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values
+          epsilon = (rmax - rmin)/(2**10)
+          error_out(ig,l) = 
+     &    100*(epsilon**3 +3*epsilon**2*rdicho +3*epsilon*rdicho**2)
+     &    / (rdicho**3-rccn**3)     
+        endif 
+        print*, 'error masse glace %', error_out(ig,l)
+        print*, 'rice,ice,vap bf',
+     &   rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap)
+       endif
+          
+c         If the initial condition is subsaturated and there is not enough ice available for sublimation 
+c         to reach equilibrium, req is neagtive. Therefore, enforce physical value.
+          rice(ig,l) = max(rdicho,rccn)
+          
+!!!!! Water ice mass is computed with radius at t+1 and their current number 
+!!!!! Nccn is at t or t+1, depending on what has been done before
+!          zq(ig,l,igcm_h2o_ice) = 
+!     &       (pi*rho_ice*zq(ig,l,igcm_ccn_number)*4./3. 
+!     &         * rdicho*rdicho*rdicho - 
+!     &         zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust)
+!     &       * tauscaling(ig)
+
+!!!!! Water ice mass is computed with radius at t+1 and their number at t+1
+!!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before
+!!!!! TO BE CLEANED ONE DAY
+          var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
+          var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig)   + dM
+          zq(ig,l,igcm_h2o_ice) = 
+     &       (pi*rho_ice*var1*4./3. 
+     &         * rdicho*rdicho*rdicho - 
+     &         var2*rho_ice/rho_dust)
+      
+     
+      !!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1))
+          zq(ig,l,igcm_h2o_ice) = 
+     &     min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap), 
+     &     zq(ig,l,igcm_h2o_ice))
+          zq(ig,l,igcm_h2o_ice) = 
+     &     max(1e-30,zq(ig,l,igcm_h2o_ice))
+     
+          zq(ig,l,igcm_h2o_vap) = 
+     &         zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
+     &       - zq(ig,l,igcm_h2o_ice)
+     
+     
+       if (test_flag) then
+        print*, 'rice,ice,vap af',
+     &     rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap)
+        print*, 'satu bf, af',
+     &     zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l),
+     &     zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
+       endif
+     
+     
+!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
+          if ((zq(ig,l,igcm_h2o_ice).le. -1e-8)
+     &    .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then
+           print*, 'NEGATIVE WATER'
+           print*, 'ig,l', ig,l
+           print*, 'satu', satu
+           print*, 'vap, ice bf', 
+     &          zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice)
+           print*, 'vap, ice af', 
+     &          zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice)
+           print*, 'ccn N,M bf',
+     &          zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass)
+           print*, 'ccn N,M af',
+     &          zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass)
+           print*, 'tauscaling',
+     &          tauscaling(ig) 
+           print*, 'req,rccn,rice bf,rice af',
+     &         req,rccn,var1,rice(ig,l)
+           print*, 'rmin, rmax', var2,var3
+           print*, 'error_out,tdicho,ptimestep',
+     &          error_out(ig,l),tdicho,ptimestep
+           print*, ' '
+          endif
+!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
+          
+          
+        ENDIF !of if Nccn >1
+          
+     
+!=============================================================
+! 5. Dust cores released, tendancies, latent heat, etc ...
+!=============================================================
+
 c         If all the ice particles sublimate, all the condensation
-c           nuclei are released:
-          if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
-c           for coherence 
-            dM = 0
-            dN = 0
-            dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
-            dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
-c           Water ice particles
+c         nuclei are released:
+          if (zq(ig,l,igcm_h2o_ice).le.1e-10) then
+!           for coherence 
+!            dM = 0
+!            dN = 0
+!            dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
+!            dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
+c           Water
+            zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) 
+     &                            + zq(ig,l,igcm_h2o_ice)
             zq(ig,l,igcm_h2o_ice) = 0.
 c           Dust particles
-            zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + 
-     &        zq(ig,l,igcm_ccn_mass)
-            zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + 
-     &        zq(ig,l,igcm_ccn_number)
+            zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
+     &                              + zq(ig,l,igcm_ccn_mass)
+            zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
+     &                                + zq(ig,l,igcm_ccn_number)
 c           CCNs
             zq(ig,l,igcm_ccn_mass) = 0.
             zq(ig,l,igcm_ccn_number) = 0.
           endif
-
-        dN = dN/ tauscaling(ig)
-        dM = dM/ tauscaling(ig)
-c       Dust particles
-        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
-        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
-c       CCNs
-        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
-        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
+          
+
+!        dN = dN/ tauscaling(ig)
+!        dM = dM/ tauscaling(ig)
+!c       Dust particles
+!        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
+!        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
+!c       CCNs
+!        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
+!        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
 
                
@@ -456,4 +622,6 @@
      
         count = count +1
+        
+        
         ELSE ! for coherence (rdust, rice computations etc ..)
           zq(ig,l,igcm_dust_mass)   = zq0(ig,l,igcm_dust_mass)
@@ -463,19 +631,23 @@
           zq(ig,l,igcm_h2o_ice)     = zq0(ig,l,igcm_h2o_ice)
           zq(ig,l,igcm_h2o_vap)     = zq0(ig,l,igcm_h2o_vap)
-
+          
 ! pour les sorties de test
-          satu_out(ig,l) = satu
-          Mcon_out(ig,l) = 0
-          newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
-          gr_out(ig,l) = 0
-          dN_out(ig,l) = 0
-          dM_out(ig,l) = 0
+!          satu_out(ig,l) = satu
+!          gr_out(ig,l) = 0
+!          dN_out(ig,l) = 0
+!          dM_out(ig,l) = 0
          
         ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
         
-c-----update temperature        
+!        ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig)
+!        ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig)
+!
+!        satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l)
+!        satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
+        
+        
+c-----update temperature (latent heat relase)     
           lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
-          pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
-c          pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
+          pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
           
 c----- update rice & rhocloud AFTER microphysic
@@ -488,5 +660,4 @@
           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
           
-          rice_out(ig,l)=rice(ig,l)
           if ((Mo.lt.1.e-20) .or. (No.le.1)) then
               rice(ig,l) = 1.e-8
@@ -495,5 +666,4 @@
      &          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)
           
           nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
@@ -510,102 +680,63 @@
           rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
 
-        ENDDO
-      ENDDO
-      
-!      print*, 'SATU'
-!      print*, satu_out(1,:)
-     
-c------------------------------------------------------------------
-c------------------------------------------------------------------
-c------------------------------------------------------------------
-c------------------------------------------------------------------
-c------------------------------------------------------------------
-c     TESTS
-
-      IF (output_sca) then
- 
+        ENDDO ! of ig loop
+      ENDDO ! of nlayer loop
+      
+     
+!!!!!!!!!!!!!! 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
+         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
+       ENDDO
+       ENDDO
+
       print*, 'count is ',count, ' i.e. ',
      &     count*100/(nlay*ngrid), '% for microphys computation'      
 
-      dM_col(:)    = 0
-      dN_col(:)    = 0
-      Mdust_col(:) = 0
-      Ndust_col(:) = 0
-      Mccn_col(:)  = 0
-      Nccn_col(:)  = 0
-      dMice_col(:) = 0
-      drice_col(:) = 0
-      icetot(:)    = 0
-      do l=1, nlay
-        do ig=1,ngrid
-          dM_col(ig) = dM_col(ig) + 
-     &       dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
-          dN_col(ig) = dN_col(ig) + 
-     &       dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
-          Mdust_col(ig) = Mdust_col(ig) + 
-     &       zq(ig,l,igcm_dust_mass)*tauscaling(ig)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          Ndust_col(ig) = Ndust_col(ig) + 
-     &       zq(ig,l,igcm_dust_number)*tauscaling(ig)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          Mccn_col(ig) = Mccn_col(ig) + 
-     &       zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          Nccn_col(ig) = Nccn_col(ig) + 
-     &       zq(ig,l,igcm_ccn_number)*tauscaling(ig)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          dMice_col(ig) = dMice_col(ig) + 
-     &       Mcon_out(ig,l)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          drice_col(ig) = drice_col(ig) + 
-     &       rice_out(ig,l)*zq(ig,l,igcm_h2o_ice)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-          icetot(ig) = icetot(ig) + 
-     &       zq(ig,l,igcm_h2o_ice)
-     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
-        enddo ! of do ig=1,ngrid
-      enddo ! of do l=1,nlay
-      
-      drice_col(ig) = drice_col(ig)/icetot(ig)
-
       IF (ngrid.ne.1) THEN ! 3D
-         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
-     &                    satu_out)
-         call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
-     &                    dM_col)
-         call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
-     &                    dN_col)
-         call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
-     &                    Ndust_col)
-         call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
-     &                    Mdust_col)
-         call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
-     &                    Nccn_col)
-         call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
-     &                    Mccn_col)
-         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
-     &                    dM_out)
-         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
-     &                    dN_out)
+!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
+!     &                    satu_out)
+!         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
+!     &                    dM_out)
+!         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
+!     &                    dN_out)
+         call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
+     &                    error2d)
+!         call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
+!     &                    zqsat)
       ENDIF
 
       IF (ngrid.eq.1) THEN ! 1D
-
-         call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
-     &                    newvap_out)
-         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
-     &                    gr_out)
-         call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
-     &                    rate_out)
-         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
-     &                    dM_out)
-         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
-     &                    dN_out)
-         call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0,
-     &                    dMice_col)
-         call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0,
-     &                    drice_col)
-         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
-     &                    Mcon_out)
+         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)
+c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
+c     &                    gr_out)
+c         call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
+c     &                    rate_out)
+c         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
+c     &                    dM_out)
+c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
+c     &                    dN_out)
          call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
      &                    zqsat)
@@ -620,19 +751,55 @@
          call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
      &                    rhocloud)
-         call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
-     &                    dM_col)
-         call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
-     &                    dN_col)
-         call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
-     &                    Ndust_col)
-         call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
-     &                    Mdust_col)
-         call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
-     &                    Nccn_col)
-         call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
-     &                    Mccn_col)
       ENDIF
-      ENDIF ! endif output_sca
-c------------------------------------------------------------------
+      
+      ENDIF ! endif test_flag
+!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
+!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
+!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
+
       return
       end
+      
+      
+      
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
+c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
+c It is an analytical solution to the ice radius growth equation, 
+c with the approximation of a constant 'reduced' cunningham correction factor 
+c (lambda in growthrate.F) taken at radius req instead of rice    
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      subroutine phi(rice,req,coeff1,coeff2,time)
+      
+      implicit none
+      
+      ! inputs
+      real rice ! ice radius
+      real req  ! ice radius at equilibirum
+      real coeff1  ! coeff for the log
+      real coeff2  ! coeff for the arctan
+
+      ! output      
+      real time
+      
+      !local
+      real var
+      
+      ! 1.73205 is sqrt(3)
+      
+      var = max(
+     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
+            
+       time = 
+     &   coeff1 * 
+     &   log( var )
+     & + coeff2 * 1.73205 *
+     &   atan( (2*rice+req) / (1.73205*req) )
+      
+      return
+      end
+      
+      
+      
Index: trunk/LMDZ.MARS/libf/phymars/lwu.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/lwu.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/lwu.F	(revision 626)
@@ -27,4 +27,6 @@
 c
 c MODIF : FF : removing the monster bug on water ice clouds 11/2010
+c
+c MODIF : TN : corrected bug if very big water ice clouds 04/2012
 c-----------------------------------------------------------------------
  
@@ -201,10 +203,14 @@
 c   et pourquoi pas d'abord?  hourdin@lmd.ens.fr
 
-           aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl))
+c TN 04/12 : if very big water ice clouds, aer_t is strictly rounded 
+c to zero in lower levels, which is a source of NaN
+           !aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl))
+           aer_t(jl,ja,jkl)=max(exp(-pview*aer_a(jl,ja,jkl)),1e-30) 
+           
           
           enddo
         enddo
-      enddo                 
-
+      enddo      
+      
 c----------------------------------------------------------------------
       return
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 626)
@@ -194,6 +194,4 @@
       REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
       
-      REAL watercapflag(ngridmx)     ! water cap flag
-
 c     Variables used by the water ice microphysical scheme:
       REAL rice(ngridmx,nlayermx)    ! Water ice geometric mean radius (m)
@@ -986,5 +984,5 @@
            call watercloud(ngrid,nlayer,ptimestep,
      &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
-     &                pq,pdq,zdqcloud,zdqscloud,zdtcloud,
+     &                pq,pdq,zdqcloud,zdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud)
@@ -1008,10 +1006,14 @@
      &                    pdq(ig,l,igcm_h2o_ice)+
      &                    zdqcloud(ig,l,igcm_h2o_ice)
-                 if (((pq(ig,l,igcm_h2o_ice) +
+             ! There are negative values issues with some tracers (17/04)
+                 if (((pq(ig,l,igcm_h2o_ice) + 
      &                 pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0)
-     &           then
-                   pdq(ig,l,igcm_h2o_ice) = 
-     &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30
-                 endif
+     &             pdq(ig,l,igcm_h2o_ice) =                
+     &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30  
+                 if (((pq(ig,l,igcm_h2o_vap) + 
+     &                 pdq(ig,l,igcm_h2o_vap)*ptimestep)) .le. 0)
+     &             pdq(ig,l,igcm_h2o_vap) =                
+     &                 - pq(ig,l,igcm_h2o_vap)/ptimestep + 1e-30
+     
                  IF (scavenging) THEN
                    pdq(ig,l,igcm_ccn_mass)=
@@ -1052,10 +1054,4 @@
            ENDIF ! of IF (water) THEN
 
-! Increment water ice surface tracer tendency
-         DO ig=1,ngrid
-           dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
-     &                               zdqscloud(ig,igcm_h2o_ice)
-         ENDDO
-         
          END IF  ! of IF (water)
 
@@ -1823,13 +1819,4 @@
      &                       'surface h2o_ice',
      &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
-
-c            if (caps) then
-c             do ig=1,ngridmx
-c                if (watercaptag(ig)) watercapflag(ig) = 1
-c             enddo
-c            CALL WRITEDIAGFI(ngridmx,'watercaptag',
-c     &                         'Ice water caps',
-c     &                         '',2,watercapflag)
-cs            endif
 c           CALL WRITEDIAGFI(ngridmx,'albedo',
 c    &                         'albedo',
Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 626)
@@ -69,7 +69,7 @@
       REAL, PARAMETER :: ccn0 = 1.3E8
 
-      LOGICAL firstcall
-      DATA firstcall/.true./
-      SAVE firstcall
+c      LOGICAL firstcall
+c      DATA firstcall/.true./
+c      SAVE firstcall
 
       REAL CBRT
@@ -83,7 +83,10 @@
 c==================================================================
 
-      IF (firstcall) THEN
+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       1.1 Dust particles
@@ -120,6 +123,7 @@
           ENDDO
         ENDIF ! of if (water.AND.activice)
-        firstcall = .false.
-      ENDIF ! of if firstcall
+        
+c        firstcall = .false.
+c      ENDIF ! of if firstcall
 
 c==================================================================
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 626)
@@ -39,4 +39,5 @@
 #include "comgeomfi.h"
 #include "tracer.h"
+#include "microphys.h"
 
 c
@@ -99,4 +100,6 @@
       LOGICAL firstcall
       SAVE firstcall
+
+      REAL lw
 
 c     variable added for CO2 condensation:
@@ -796,4 +799,9 @@
 c             Starting upward calculations for water :
                zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
+c             Take into account H2o latent heat in surface energy budget          
+               lw = (2834.3-0.28*(ptsrf(ig)-To)
+     &              -0.004*(ptsrf(ig)-To)**2)*1.e+3
+               pdtsrf(ig) = pdtsrf(ig) 
+     &                    + pdqsdif(ig,igcm_h2o_ice)*lw/pcapcal(ig)
             ENDDO ! of DO ig=1,ngrid
           ELSE
@@ -806,5 +814,5 @@
           DO ilay=2,nlay
              DO ig=1,ngrid
-                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
              ENDDO
           ENDDO
Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 625)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 626)
@@ -1,9 +1,7 @@
-       SUBROUTINE watercloud(ngrid,nlay,ptimestep, 
+       SUBROUTINE watercloud(ngrid,nlay, ptimestep, 
      &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
-     &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
+     &                pq,pdq,pdqcloud,pdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud)
-! to use  'getin'
-      USE ioipsl_getincom
       IMPLICIT NONE
 
@@ -15,11 +13,8 @@
 c    - An improved microphysical scheme (see improvedclouds.F)
 c
-c  There is a time loop specific to cloud formation and sedimentation
-c  due to timescales smaller than the GCM integration timestep.
-c
 c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 
 c           J.-B. Madeleine, Thomas Navarro
 c
-c  2004 - 2012
+c  2004 - Oct. 2011
 c=======================================================================
 
@@ -40,28 +35,27 @@
 
       INTEGER ngrid,nlay
-      INTEGER nq                 ! nombre de traceurs 
+      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 pdpsrf(ngrid)         ! tendence surf pressure
+      REAL pdpsrf(ngrid)         ! tendance surf pressure
       REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
       REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
       REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
-      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
+      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
 
       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
-      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
+      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
 
-      REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
-      REAL tauscaling(ngridmx)   ! Convertion factor for dust amount
-      real rdust(ngridmx,nlay)   ! Dust geometric mean radius (m)
+      REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
+      REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
+      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
 
 c   Outputs:
 c   -------
 
-      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
-      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
-      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
-                                   ! a la chaleur latente
+      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
+      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
+                                   !   a la chaleur latente
 
       REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
@@ -69,42 +63,11 @@
       REAL nuice(ngrid,nlay)   ! Estimated effective variance
                                !   of the size distribution
-      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
-      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
+      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
+      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
 
 c   local:
 c   ------
 
-      ! for sedimentation
-      REAL zq(ngridmx,nlay,nqmx)    ! local value of tracers
-      real masse (ngridmx,nlay)     ! Layer mass (kg.m-2)
-      real epaisseur (ngridmx,nlay) ! Layer thickness (m)
-      real wq(ngridmx,nlay+1)       ! displaced tracer mass (kg.m-2)
-      
-      ! for ice radius computation
-      REAL Mo,No
-      REAl ccntyp
-      real beta      ! correction for the shape of the ice particles (cf. newsedim)
-      save beta
-      
-      ! for time loop
-      INTEGER microstep  ! time subsampling step variable
-      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
-      SAVE imicro
-      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
-      SAVE microtimestep
-      
-      ! tendency given by clouds (inside the micro loop)
-      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
-      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
-
-      ! global tendency (clouds+sedim+physics)
-      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
-      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
-      
-      REAL CBRT
-      EXTERNAL CBRT
-
-      
-      INTEGER iq,ig,l
+      INTEGER ig,l
       LOGICAL,SAVE :: firstcall=.true.
 
@@ -129,301 +92,22 @@
         write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
         write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
-                
-
-        write(*,*) "correction for the shape of the ice particles ?"
-        beta=0.75 ! default value
-        call getin("ice_shape",beta)
-        write(*,*) " ice_shape = ",beta
-        
-        write(*,*) "time subsampling for microphysic ?"
-        imicro = 1
-        call getin("imicro",imicro)
-        write(*,*)"imicro = ",imicro
-        
-        microtimestep = ptimestep/float(imicro)
-        write(*,*)"Physical timestep is",ptimestep 
-        write(*,*)"Microphysics timestep is",microtimestep 
 
         firstcall=.false.
       ENDIF ! of IF (firstcall)
       
-c-----Initialization
-      subpdq(1:ngrid,1:nlay,1:nq) = 0
-      subpdt(1:ngrid,1:nlay)      = 0
-      pdqscloud(1:ngrid,1:nq)     = 0
-      zq(1:ngrid,1:nlay,1:nq)     = 0
 
-c-----Computing the different layer properties for clouds sedimentation      
-      do l=1,nlay
-        do ig=1, ngrid
-          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
-          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
-         enddo
-      enddo
-
-c------------------------------------------------------------------
-c Time subsampling for coupled microphysics and sedimentation
-c------------------------------------------------------------------
-      DO microstep=1,imicro 
-      
-c-------------------------------------------------------------------
-c   1.  Tendencies: 
-c------------------
-
-
-c------ Temperature tendency subpdt
-        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
-        ! If imicro=1 subpdt is the same as pdt
-        DO l=1,nlay
-          DO ig=1,ngrid
-             subpdt(ig,l) = subpdt(ig,l)
-     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
-          ENDDO
-        ENDDO
-c------ Traceurs tendencies subpdq
-c------ At each micro timestep we add pdq in order to have a stepped entry
-        IF (microphys) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              subpdq(ig,l,igcm_dust_mass) = 
-     &            subpdq(ig,l,igcm_dust_mass)
-     &          + pdq(ig,l,igcm_dust_mass)
-              subpdq(ig,l,igcm_dust_number) = 
-     &            subpdq(ig,l,igcm_dust_number)
-     &          + pdq(ig,l,igcm_dust_number)
-              subpdq(ig,l,igcm_ccn_mass) = 
-     &            subpdq(ig,l,igcm_ccn_mass)
-     &          + pdq(ig,l,igcm_ccn_mass)
-              subpdq(ig,l,igcm_ccn_number) = 
-     &            subpdq(ig,l,igcm_ccn_number)
-     &          + pdq(ig,l,igcm_ccn_number)
-            ENDDO
-          ENDDO
-        ENDIF
-        DO l=1,nlay
-          DO ig=1,ngrid
-            subpdq(ig,l,igcm_h2o_ice) = 
-     &          subpdq(ig,l,igcm_h2o_ice)
-     &        + pdq(ig,l,igcm_h2o_ice)
-            subpdq(ig,l,igcm_h2o_vap) = 
-     &          subpdq(ig,l,igcm_h2o_vap)
-     &        + pdq(ig,l,igcm_h2o_vap)
-          ENDDO
-        ENDDO
-        
-        
-c-------------------------------------------------------------------
-c   2.  Main call to the different cloud schemes:
-c------------------------------------------------
-        IF (microphys) THEN
-           CALL improvedclouds(ngrid,nlay,microtimestep,
-     &             pplev,pplay,pzlev,pt,subpdt,
-     &             pq,subpdq,subpdqcloud,subpdtcloud,
+c     Main call to the different cloud schemes:
+      IF (microphys) THEN
+        CALL improvedclouds(ngrid,nlay,ptimestep,
+     &             pplev,pplay,pt,pdt,
+     &             pq,pdq,pdqcloud,pdtcloud,
      &             nq,tauscaling,rdust,rice,nuice,
      &             rsedcloud,rhocloud)
-        ELSE
-           CALL simpleclouds(ngrid,nlay,microtimestep,
-     &             pplev,pplay,pzlev,pzlay,pt,subpdt,
-     &             pq,subpdq,subpdqcloud,subpdtcloud,
+      ELSE
+        CALL simpleclouds(ngrid,nlay,ptimestep,
+     &             pplev,pplay,pzlev,pzlay,pt,pdt,
+     &             pq,pdq,pdqcloud,pdtcloud,
      &             nq,tau,rice,nuice,rsedcloud)
-        ENDIF
-      
-c--------------------------------------------------------------------
-c    3.  CCN & ice sedimentation:
-c--------------------------------
-! Sedimentation is done here for water ice and its CCN only
-! callsedim in physics is done for dust (and others species if any)           
-                   
-        DO l=1,nlay
-         DO ig=1,ngrid
-          subpdt(ig,l) =
-     &      subpdt(ig,l) + subpdtcloud(ig,l)
-         ENDDO
-        ENDDO
-        
-c------- water ice update before sedimentation (radius is done in the cloud scheme routine)    
-         DO l=1,nlay
-          DO ig=1, ngrid
-          zq(ig,l,igcm_h2o_ice)= max(1e-30,
-     &       pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice)
-     &       + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep)
-!          zq(ig,l,igcm_h2o_vap)= max(1e-30,
-!     &       pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap)
-!     &       + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep)
-          ENDDO
-         ENDDO
-        
-      
-c------- CCN update before sedimentation  
-        IF (microphys) THEN
-         DO l=1,nlay
-          DO ig=1,ngrid
-           zq(ig,l,igcm_ccn_number)=
-     &       pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number)
-     &       + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
-           zq(ig,l,igcm_ccn_number)=  max(1e-30,
-     &       zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ?
-           zq(ig,l,igcm_ccn_mass)=
-     &       pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass)
-     &       + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
-           zq(ig,l,igcm_ccn_mass)=max(1e-30,
-     &       zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ?
-          ENDDO
-         ENDDO
-        ENDIF
-        
-
-        
-        IF (microphys) THEN
-        
-c------- CCN number sedimentation
-          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
-     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
-     &        rhocloud,zq(1,1,igcm_ccn_number),wq,beta)
-          do ig=1,ngrid
-            ! matters if one would like to know ccn surface deposition
-            pdqscloud(ig,igcm_ccn_number)=
-     &          pdqscloud(ig,igcm_ccn_number) + wq(ig,1)
-          enddo
-          
-c------- CCN mass sedimentation
-          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
-     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
-     &        rhocloud,zq(1,1,igcm_ccn_mass),wq,beta)
-          do ig=1,ngrid
-            ! matters if one would like to know ccn surface deposition
-            pdqscloud(ig,igcm_ccn_mass)=
-     &          pdqscloud(ig,igcm_ccn_mass) + wq(ig,1)
-          enddo
-          
-c------- Water ice sedimentation -- improved microphys
-          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
-     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
-     &        rhocloud,zq(1,1,igcm_h2o_ice),wq,beta)
-        ELSE
-        
-c------- Water ice sedimentation -- simple microphys
-          call newsedim(ngrid,nlay,ngrid*nlay,1,
-     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
-     &        rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta)
-     
-        ENDIF
-
-
-c------- Surface ice tendency update
-        DO ig=1,ngrid
-          pdqscloud(ig,igcm_h2o_ice)=
-     &          pdqscloud(ig,igcm_h2o_ice) + wq(ig,1)
-        ENDDO
-      
-
-c-------------------------------------------------------------------
-c   5.  Updating tendencies after sedimentation:
-c-----------------------------------------------
-
-        DO l=1,nlay
-         DO ig=1,ngrid
-           
-           subpdq(ig,l,igcm_h2o_ice) = 
-     &     (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice))
-     &      /microtimestep
-
-     
-           subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap)
-     &        +subpdqcloud(ig,l,igcm_h2o_vap)
-     
-         ENDDO
-        ENDDO
-      
-      
-        IF (microphys) then
-         DO l=1,nlay
-          DO ig=1,ngrid
-           subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
-     &        -pq(ig,l,igcm_ccn_number))/microtimestep
-           subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
-     &        -pq(ig,l,igcm_ccn_mass))/microtimestep
-          ENDDO
-         ENDDO
-        ENDIF
-        
-
-
-      ENDDO ! of DO microstep=1,imicro
-      
-c-------------------------------------------------------------------
-c   6.  Compute final tendencies after time loop:
-c------------------------------------------------
-
-c------ Whole temperature tendency pdtcloud
-       DO l=1,nlay
-         DO ig=1,ngrid
-             pdtcloud(ig,l) =
-     &         subpdt(ig,l)/imicro-pdt(ig,l)
-          ENDDO
-       ENDDO
-c------ Traceurs tendencies pdqcloud
-       DO iq=1,nq
-        DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro
-     &       - pdq(ig,l,iq)
-         ENDDO
-        ENDDO
-       ENDDO
-c------ Traceurs surface tendencies pdqscloud
-       DO iq=1,nq
-        DO ig=1,ngrid
-           pdqscloud(ig,iq) =
-     &         pdqscloud(ig,iq)/ptimestep
-        ENDDO
-       ENDDO
-       
-
-          
-c------Update the ice particle size "rice" for output or photochemistry.
-c------It is not used again in the water cycle.
-       IF(scavenging) THEN 
-        DO l=1, nlay
-         DO ig=1,ngrid
-          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.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
-       ELSE
-        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))
-          ENDDO
-        ENDDO
-       ENDIF ! of IF(scavenging)
-                
-!--------------------------------------------------------------
-!--------------------------------------------------------------
+      ENDIF
       
 
@@ -439,6 +123,5 @@
       end do
 
-c=======================================================================
-
+      RETURN
       END
  
