Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 632)
+++ /trunk/LMDZ.MARS/README	(revision 633)
@@ -1622,2 +1622,5 @@
 >> Some cleanup on messages and comments in code about the reference pressure
    for dust opacity which is now 610 Pa.
+
+== 25/04/12 == TN
+>> The new scheme does not work for now. Instead, use of a subtimestep for nucleation and growth.
Index: /trunk/LMDZ.MARS/libf/phymars/growthrate.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 633)
@@ -1,3 +1,3 @@
-      subroutine growthrate(temp,pmid,psat,seq,rcrystal,coeff1,coeff2)
+      subroutine growthrate(temp,pmid,psat,rcrystal,res)
 
       IMPLICIT NONE
@@ -31,10 +31,9 @@
       REAL temp     ! temperature in the middle of the layer (K)
       REAL pmid     ! pressure in the middle of the layer (K)
-      REAL*8 psat   ! water vapor saturation pressure (Pa) 
-      REAL*8 seq    ! Equilibrium saturation ratio
+      REAL psat   ! water vapor saturation pressure (Pa) 
       REAL rcrystal ! crystal radius before condensation (m)
 
 c     Output
-      REAL coeff1,coeff2     ! coefficients for the analytical relation between time and radius
+      REAL res      ! growth resistance (res=Rk+Rd)
 
 
@@ -42,8 +41,8 @@
 c   ------
 
-      REAL*8 k,Lv                 
-      REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
-      REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
-      REAL*8 Rk,Rd
+      REAL k,Lv                 
+      REAL knudsen           ! Knudsen number (gas mean free path/particle radius)
+      REAL afactor,Dv,lambda       ! Intermediate computations for growth rate
+      REAL Rk,Rd
       
       
@@ -86,11 +85,13 @@
 
 c     - Compute Rk
-      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp*temp)
+      Rk = Lv*Lv* rho_ice * mh2o / (k*rgp*temp*temp)
 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))
+      res = Rk + Rd*(1. + lambda * knudsen)
+      
+      !coeff1 = real(Rk + Rd*(1. + lambda * knudsen))
+      !coeff2 = real(Rk + Rd*(1. - lambda * knudsen))
       
 c Below are growth rate used for other schemes
Index: /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 633)
@@ -1,10 +1,11 @@
       subroutine improvedclouds(ngrid,nlay,ptimestep,
-     &             pplev,pplay,pt,pdt,
+     &             pplay,pt,pdt,
      &             pq,pdq,pdqcloud,pdtcloud,
-     &             nq,tauscaling,rdust,rice,nuice,
-     &             rsedcloud,rhocloud)
+     &             nq,tauscaling)
 ! to use  'getin'
       USE ioipsl_getincom
       implicit none
+      
+      
 c------------------------------------------------------------------
 c  This routine is used to form clouds when a parcel of the GCM is
@@ -21,6 +22,4 @@
 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
@@ -43,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)
             
@@ -54,13 +52,6 @@
                                  !   (kg/kg.s-1)
       REAL tauscaling(ngridmx)     ! Convertion factor for qdust and Ndust
-      REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
 
 c     Outputs:
-      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
-                                 ! (r_c in montmessin_2004)
-      REAL nuice(ngrid,nlay)     ! Estimated effective variance
-                                 !   of the size distribution
-      REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
-      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
       REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
                                    !   H2O(kg/kg.s-1)
@@ -89,19 +80,28 @@
       REAL zqsat(ngridmx,nlayermx)    ! saturation
       REAL lw                         !Latent heat of sublimation (J.kg-1) 
-      REAL Cste
-      REAL dMice           ! mass of condensated ice
-      REAL sumcheck
+      REAL cste
+      REAL dMice           ! mass of condensed ice
+!      REAL sumcheck
       REAL*8 ph2o          ! Water vapor partial pressure (Pa)
       REAL*8 satu          ! Water vapor saturation ratio over ice
       REAL*8 Mo,No
-      REAL*8 dN,dM,newvap
-      REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
+      REAL*8 Rn, Rm, dev2, 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
-      REAl*8 seq
+
       REAL*8 sig      ! Water-ice/air surface tension  (N.m)
       EXTERNAL sig
+
+      REAL dN,dM
+      REAL rate(nbin_cld)  ! nucleation rate
+      REAL seq
+
+      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
+                                 ! (r_c in montmessin_2004)
+      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
+      REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
+
+      REAL res      ! Resistance growth
+
 
 c     Parameters of the size discretization
@@ -118,34 +118,18 @@
       DOUBLE PRECISION vol_cld(nbin_cld)  ! particle volume for each bin (m3)
 
+
       REAL sigma_ice ! Variance of the ice and CCN distributions
       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 -- 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 Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
-      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
+c TESTS
+
+      INTEGER countcells
       
       LOGICAL test_flag    ! flag for test/debuging outputs
-      
-      test_flag = .false.
-c----------------------------------      
-c----------------------------------      
+      SAVE    test_flag     
+ 
 
 c------------------------------------------------------------------
@@ -164,4 +148,6 @@
 
 c       Volume ratio between two adjacent bins
+   !     vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
+   !     vrat_cld = exp(vrat_cld)
         vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
         vrat_cld = dexp(vrat_cld)
@@ -171,4 +157,5 @@
         rad_cld(1) = rmin_cld
         vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
+   !     vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
 
         do i=1,nbin_cld-1
@@ -197,4 +184,5 @@
 
         do i=1,nbin_cld+1
+    !        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
@@ -210,74 +198,58 @@
 c       Variance of the ice and CCN distributions
         sigma_ice = sqrt(log(1.+nuice_sed))
-        
+       
+ 
         write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
         write(*,*) 'nuice for sedimentation:', nuice_sed
         write(*,*) 'Volume of a water molecule:', vo1
 
+
+        test_flag = .false.
+ 
         firstcall=.false.
       END IF
 
+
 !=============================================================
 ! 1. Initialisation
 !=============================================================
+      cste = 4*pi*rho_ice*ptimestep
 
 c     Initialize the tendencies
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
       pdtcloud(1:ngrid,1:nlay)=0
-      
-c     Update the needed variables
-      do l=1,nlay
-        do ig=1,ngrid
-c         Temperature
-          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
-c         Dust mass mixing ratio
-          zq(ig,l,igcm_dust_mass) = 
-     &      pq(ig,l,igcm_dust_mass) + 
-     &      pdq(ig,l,igcm_dust_mass) * ptimestep
-          zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
-c         Dust particle number
-          zq(ig,l,igcm_dust_number) = 
-     &      pq(ig,l,igcm_dust_number) + 
-     &      pdq(ig,l,igcm_dust_number) * ptimestep
-          zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
-c         Update rdust from last tendencies
-          rdust(ig,l)=
-     &      CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
-     &      max(zq(ig,l,igcm_dust_number),0.01))
-          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)         
-c         CCN mass mixing ratio
-          zq(ig,l,igcm_ccn_mass)=
-     &      pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
-          zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
-          zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
-c         CCN particle number
-          zq(ig,l,igcm_ccn_number)=
-     &      pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
-          zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
-          zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
-c         Water vapor
-          zq(ig,l,igcm_h2o_vap)=
-     &      pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
-          zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004 
-          zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
-c         Water ice
-          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),1e-30) ! FF 12/2004 
-          zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
-        enddo
-      enddo
-
+    
+c     Initialize the tendencies
+      pdqcloud(1:ngrid,1:nlay,1:nq)=0
+      pdtcloud(1:ngrid,1:nlay)=0
+      
+c     Initialize the tendencies
+      pdqcloud(1:ngrid,1:nlay,1:nq)=0
+      pdtcloud(1:ngrid,1:nlay)=0
+      
+      
+      zt(1:ngrid,1:nlay) = 
+     &      pt(1:ngrid,1:nlay) + 
+     &      pdt(1:ngrid,1:nlay) * ptimestep
+
+      zq(1:ngrid,1:nlay,1:nq) = 
+     &      pq(1:ngrid,1:nlay,1:nq) + 
+     &      pdq(1:ngrid,1:nlay,1:nq) * ptimestep
+      
+      zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
+      
+      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
+     &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30
+      
 !=============================================================
 ! 2. Compute saturation
 !=============================================================
 
+
       dev2 = 1. / ( sqrt(2.) * sigma_ice )
-      error_out(:,:) = 0.
-
 
       call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
             
-      count = 0
+      countcells = 0
 
 c     Main loop over the GCM's grid
@@ -288,8 +260,4 @@
         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. )       .or.                             ! if there is condensation
-     &      ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
-
 
 !=============================================================
@@ -297,6 +265,15 @@
 !=============================================================
 
+       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)
+
+
 c       Expand the dust moments into a binned distribution
-        Mo = zq(ig,l,igcm_dust_mass)* 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)
@@ -313,16 +290,4 @@
           m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
         enddo
-!!! MORE EFFICIENT COMPUTATIONALLY THAN
-!        Rn = rdust(ig,l)
-!        Rm = Rn * exp( 3. * sigma_ice**2. )
-!        Rn = 1. / Rn
-!        Rm = 1. / Rm
-!        do i = 1, nbin_cld
-!          n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 )
-!     &                           -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
-!          m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 )
-!     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
-!        enddo
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         
 !        sumcheck = 0
@@ -349,27 +314,25 @@
 !          print*, "Dust binned distribution", m_aer
 !        endif
-                
+
+  
 c       Get the rates of nucleation
         call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
         
-
         dN = 0.
         dM = 0.
-        rate_out(ig,l) = 0
         do i = 1, nbin_cld
-          n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
-          m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
+          n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
+          m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
           dN       = dN + n_aer(i) * rate(i) * ptimestep
           dM       = dM + m_aer(i) * rate(i) * ptimestep
-          !rate_out(ig,l)=rate_out(ig,l)+rate(i)
         enddo
 
-c       Update CCNs, can also be done after the radius growth
-c       Dust particles
+
+c       Update 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
+c       Update CCNs
         zq(ig,l,igcm_ccn_mass)   = 
      &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
@@ -377,4 +340,5 @@
      &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
         
+        ENDIF ! of is satu >1
 
 !=============================================================
@@ -382,7 +346,15 @@
 !=============================================================
 
+c We trigger crystal growth if and only if there is at least one nuclei (N>1).
+c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
+c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
+
+       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
-        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
+     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
+        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
+
 
         rhocloud(ig,l) =  zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
@@ -391,14 +363,9 @@
         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)   
+     &      CBRT( real(Mo/No) * 0.75 / pi / rhocloud(ig,l) )
+
 
 c       saturation at equilibrium
@@ -406,167 +373,33 @@
      &           (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
-            !  crystal size is decreasing
-            rmax = max(min(rice(ig,l),req),rccn) 
-            rmin = max(rice(ig,l),req)
-          endif
-                         !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           
+c       get resistance growth
+        call growthrate(zt(ig,l),pplay(ig,l),
+     &             real(ph2o/satu),rice(ig,l),res)
+
+
+ccccccc  implicit scheme of mass growth
+
+        dMice =
+     & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
+     &   /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
+
+
+! With the above scheme, dMice cannot be bigger than vapor, 
+! but can be bigger than all available ice.
+       dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
+       dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
+
+       zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
+       zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
+
+
+       countcells = countcells + 1 
+       
+       ! latent heat release
+       lw=(2834.3-0.28*(zt(ig,l)-To)-
+     &     0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
+       pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
           
           
-        ENDIF !of if Nccn >1
-          
      
 !=============================================================
@@ -576,10 +409,6 @@
 c         If all the ice particles sublimate, all the condensation
 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)
+          if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
+          
 c           Water
             zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) 
@@ -594,92 +423,17 @@
             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
-
-               
-        pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
-     &                         -zq0(ig,l,igcm_dust_mass))/ptimestep
-        pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
-     &                         -zq0(ig,l,igcm_dust_number))/ptimestep
-        pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
-     &                         -zq0(ig,l,igcm_ccn_mass))/ptimestep
-        pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
-     &                         -zq0(ig,l,igcm_ccn_number))/ptimestep
-        pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
-     &                         -zq0(ig,l,igcm_h2o_vap))/ptimestep
-        pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
-     &                         -zq0(ig,l,igcm_h2o_ice))/ptimestep
-     
-     
-        count = count +1
-        
-        
-        ELSE ! for coherence (rdust, rice computations etc ..)
-          zq(ig,l,igcm_dust_mass)   = zq0(ig,l,igcm_dust_mass)
-          zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
-          zq(ig,l,igcm_ccn_mass)    = zq0(ig,l,igcm_ccn_mass)
-          zq(ig,l,igcm_ccn_number)  = zq0(ig,l,igcm_ccn_number)
-          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
-!          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)
-        
-!        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----- update rice & rhocloud AFTER microphysic
-          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)
-          
-          if ((Mo.lt.1.e-20) .or. (No.le.1)) then
-              rice(ig,l) = 1.e-8
-          else
-              rice(ig,l) =
-     &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
-          endif
-          
-          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
-
-c----- update rdust and sedimentation radius
-          rdust(ig,l)=
-     &      CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
-     &      max(zq(ig,l,igcm_dust_number),0.01))
-          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
-          
-          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)
-
+          ENDIF !of if Nccn>1
+         
         ENDDO ! of ig loop
       ENDDO ! of nlayer loop
+      
+      
+      ! 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
+     
       
      
@@ -689,15 +443,15 @@
       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'      
-
-      IF (ngrid.ne.1) THEN ! 3D
+!       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 ',countcells, ' i.e. ',
+     &     countcells*100/(nlay*ngrid), '% for microphys computation'
+
+!      IF (ngrid.ne.1) THEN ! 3D
 !         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
 !     &                    satu_out)
@@ -706,29 +460,29 @@
 !         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
 !     &                    dN_out)
-         call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
-     &                    error2d)
+!         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,"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)
+!      ENDIF
+
+!      IF (ngrid.eq.1) THEN ! 1D
+!         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)
@@ -739,17 +493,17 @@
 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,"satu","ratio saturation","",1,
-     &                    satu_out)
-         call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
-     &                    rice)
-         call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
-     &                    rdust)
-         call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
-     &                    rsedcloud)
-         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
-     &                    rhocloud)
-      ENDIF
+!         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,"rdust_sca","rdust","m",1,
+!     &                    rdust)
+!         call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
+!     &                    rsedcloud)
+!         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
+!     &                    rhocloud)
+!      ENDIF
       
       ENDIF ! endif test_flag
@@ -772,34 +526,34 @@
 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
-      
-      
-      
+c      subroutine phi(rice,req,coeff1,coeff2,time)
+c      
+c      implicit none
+c      
+c      ! inputs
+c      real rice ! ice radius
+c      real req  ! ice radius at equilibirum
+c      real coeff1  ! coeff for the log
+c      real coeff2  ! coeff for the arctan
+c
+c      ! output      
+c      real time
+c      
+c      !local
+c      real var
+c      
+c      ! 1.73205 is sqrt(3)
+c      
+c      var = max(
+c     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
+c            
+c       time = 
+c     &   coeff1 * 
+c     &   log( var )
+c     & + coeff2 * 1.73205 *
+c     &   atan( (2*rice+req) / (1.73205*req) )
+c      
+c      return
+c      end
+      
+      
+      
Index: /trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 633)
@@ -252,5 +252,5 @@
          nircorr=0      !default value. this is OK below 60 km.
 #else
-         nircorr=1      !default value
+         nircorr=0      !default value
 #endif
          call getin("nircorr",nircorr)
@@ -531,13 +531,5 @@
          call getin("tituscap",tituscap)
          write(*,*) "tituscap",tituscap
-
-!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
-         write(*,*) "temp tag for water caps ?"
-         temptag=.false. 
-         call getin("temptag",temptag)
-         write(*,*) " temptag = ",temptag
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!              
+                     
 
          write(*,*) "photochemistry: include chemical species"
Index: /trunk/LMDZ.MARS/libf/phymars/nuclea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/nuclea.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/nuclea.F	(revision 633)
@@ -26,5 +26,6 @@
 
 c     Output
-      DOUBLE PRECISION nucrate(nbin_cld)
+   !   DOUBLE PRECISION nucrate(nbin_cld)
+      REAL nucrate(nbin_cld)
 
 c     Local variables
@@ -119,5 +120,5 @@
             nucrate(i) = 0.
           else
-            nucrate(i)= sqrt ( fistar /
+            nucrate(i)= real(sqrt ( fistar /
      &               (3.*pi*kbz*temp*(gstar*gstar)) )
      &                  * kbz * temp * rstar
@@ -126,5 +127,5 @@
      &                  * ( nh2o*rad_cld(i) )
      &                  / ( zefshape * nus * m0 )
-     &                  * dexp (deltaf)
+     &                  * dexp (deltaf))
           endif
 
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 633)
@@ -321,5 +321,4 @@
                                !   (direct comparison with TES)
 
-
 c Test 1d/3d scavenging
       real h2otot(ngridmx)
@@ -355,4 +354,7 @@
       REAL zu2(ngridmx)
 c=======================================================================
+
+      print*,'physiq0', nq,nqmx
+
 
 c 1. Initialisation:
@@ -851,5 +853,5 @@
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
      $ dtke_th,hfmax_th,wstar)
- 
+     
          DO l=1,nlayer
            DO ig=1,ngrid
@@ -985,74 +987,56 @@
 
            call watercloud(ngrid,nlayer,ptimestep,
-     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
+     &                pplev,pplay,pdpsrf,zzlay, pt,pdt,
      &                pq,pdq,zdqcloud,zdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud)
+     
+c Temperature variation due to latent heat release
            if (activice) then
-c Temperature variation due to latent heat release
-           DO l=1,nlayer
-             DO ig=1,ngrid
-               pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
-             ENDDO
-           ENDDO
+               pdt(1:ngrid,1:nlayer) = 
+     &          pdt(1:ngrid,1:nlayer) + 
+     &          zdtcloud(1:ngrid,1:nlayer)
            endif
-
+           
 ! increment water vapour and ice atmospheric tracers tendencies
-           IF (water) THEN
-          DO l=1,nlayer
-               DO ig=1,ngrid
-                 pdq(ig,l,igcm_h2o_vap)=
-     &                     pdq(ig,l,igcm_h2o_vap)+
-     &                    zdqcloud(ig,l,igcm_h2o_vap)
-                 pdq(ig,l,igcm_h2o_ice)=
-     &                    pdq(ig,l,igcm_h2o_ice)+
-     &                    zdqcloud(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)
-     &             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)=
-     &                      pdq(ig,l,igcm_ccn_mass)+
-     &                      zdqcloud(ig,l,igcm_ccn_mass)
-                   pdq(ig,l,igcm_ccn_number)=
-     &                      pdq(ig,l,igcm_ccn_number)+
-     &                      zdqcloud(ig,l,igcm_ccn_number)
-                   pdq(ig,l,igcm_dust_mass)=
-     &                      pdq(ig,l,igcm_dust_mass)+
-     &                      zdqcloud(ig,l,igcm_dust_mass)
-                   pdq(ig,l,igcm_dust_number)=
-     &                      pdq(ig,l,igcm_dust_number)+
-     &                      zdqcloud(ig,l,igcm_dust_number)
-!!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn & Ndust > 0
-!!!!!!!!!!!!!!!!!!!!! This is due to single precision rounding problems  
-                   if (((pq(ig,l,igcm_ccn_number) +
-     &                 pdq(ig,l,igcm_ccn_number)*ptimestep)) .le. 0)
-     &             then
-                       pdq(ig,l,igcm_ccn_number) = 
-     &                 - pq(ig,l,igcm_ccn_number)/ptimestep + 1e-30
-                       pdq(ig,l,igcm_ccn_mass) = 
-     &                 - pq(ig,l,igcm_ccn_mass)/ptimestep + 1e-30
-                   endif
-                   if (((pq(ig,l,igcm_dust_number) +
-     &                 pdq(ig,l,igcm_dust_number)*ptimestep)) .le. 0)
-     &             then
-                       pdq(ig,l,igcm_dust_number) = 
-     &                 - pq(ig,l,igcm_dust_number)/ptimestep + 1e-30
-                       pdq(ig,l,igcm_dust_mass) = 
-     &                 - pq(ig,l,igcm_dust_mass)/ptimestep + 1e-30
-                   endif
-!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!
-                 ENDIF
-               ENDDO
-             ENDDO
-           ENDIF ! of IF (water) THEN
+           pdq(1:ngrid,1:nlayer,1:nq) = 
+     &       pdq(1:ngrid,1:nlayer,1:nq) + 
+     &       zdqcloud(1:ngrid,1:nlayer,1:nq)
+           
+! We need to check that we have Nccn & Ndust > 0
+! This is due to single precision rounding problems
+         if (scavenging) then
+
+           where (pq(:,:,igcm_ccn_mass) + 
+     &      ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
+              pdq(:,:,igcm_ccn_mass) = 
+     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+              pdq(:,:,igcm_ccn_number) = 
+     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+           end where
+           where (pq(:,:,igcm_ccn_number) + 
+     &      ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
+              pdq(:,:,igcm_ccn_mass) = 
+     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+              pdq(:,:,igcm_ccn_number) = 
+     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+           end where
+           where (pq(:,:,igcm_dust_mass) + 
+     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
+              pdq(:,:,igcm_dust_mass) = 
+     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+              pdq(:,:,igcm_dust_number) = 
+     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+           end where
+           where (pq(:,:,igcm_dust_number) + 
+     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
+              pdq(:,:,igcm_dust_mass) = 
+     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+              pdq(:,:,igcm_dust_number) = 
+     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+           end where
+           
+         endif ! of if scavenging
+                      
 
          END IF  ! of IF (water)
@@ -1101,6 +1085,4 @@
      
      
-      !print*, 'h2o_ice zdqsed ds physiq', zdqsed(1,:,igcm_h2o_ice)
-
            DO iq=1, nq
              DO l=1,nlayer
@@ -2121,24 +2103,12 @@
      &                'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number))
      
-        call WRITEDIAGFI(ngridmx,'zdqcloud*dt ice','cloud ice',
-     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)*ptimestep)
-        call WRITEDIAGFI(ngridmx,'zdqcloud*dt vap','cloud vap',
-     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)*ptimestep)
-        call WRITEDIAGFI(ngridmx,'zdqdif*dt ice','dif ice',
-     &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_ice)*ptimestep)
-        call WRITEDIAGFI(ngridmx,'zdqdif*dt vap','dif vap',
-     &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_vap)*ptimestep)
-        call WRITEDIAGFI(ngridmx,'zdqdif*dt vap 0','dif vap',
-     &            'kg.m-2.s-1',0,zdqdif(1,1,igcm_h2o_vap)*ptimestep)
-
-!       if (scavenging) then
-!         call WRITEDIAGFI(ngridmx,'zdqsed_ccnq','sedimentation q',
-!     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_ccn_mass))
-!         call WRITEDIAGFI(ngridmx,'zdqsed_ccnN','sedimentation N',
-!    &                 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_ccn_number))
-!      endif
-!       call WRITEDIAGFI(ngridmx,'zdqsed_ice','sedimentation q',
-!     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_h2o_ice))
-!     
+        call WRITEDIAGFI(ngridmx,'zdqcloud_ice','cloud ice',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
+        call WRITEDIAGFI(ngridmx,'zdqcloud_vap','cloud vap',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
+        call WRITEDIAGFI(ngridmx,'zdqcloud','cloud ice',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
+     &                          +zdqcloud(1,:,igcm_h2o_vap))
+     
      
           call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
Index: /trunk/LMDZ.MARS/libf/phymars/simpleclouds.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 633)
@@ -1,6 +1,6 @@
       subroutine simpleclouds(ngrid,nlay,ptimestep,
-     &             pplev,pplay,pzlev,pzlay,pt,pdt,
+     &             pplev,pplay,pzlay,pt,pdt,
      &             pq,pdq,pdqcloud,pdtcloud,
-     &             nq,tau,rice,nuice,rsedcloud)
+     &             nq,tau,rice)
       implicit none
 c------------------------------------------------------------------
@@ -43,5 +43,4 @@
       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      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
@@ -57,7 +56,4 @@
       REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
                                  ! (r_c in montmessin_2004)
-      REAL nuice(ngrid,nlay)     ! Estimated effective variance
-                                 !   of the size distribution
-      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
       real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
                                    !   H2O(kg/kg.s-1)
@@ -181,13 +177,5 @@
      &      +ccntyp(ig,l)*(4./3.)*pi*rdusttyp(ig,l)**3.)
      &      /(ccntyp(ig,l)*4./3.*pi) ), rdusttyp(ig,l)) 
-c         Effective variance of the size distribution
-          nuice(ig,l)=nuice_ref
-
-c         Sedimentation radius:
-c         ~~~~~~~~~~~~~~~~~~~~
-
-          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
-     &                         rdusttyp(ig,l) )
-          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
+
 
         enddo ! of do ig=1,ngrid
Index: /trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 633)
@@ -23,4 +23,7 @@
       REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
       REAL icedryness ! ice dryness
+      
+      ! longwatercaptag is watercaptag. Trick for some compilers
+      LOGICAL, DIMENSION(100000) :: longwatercaptag
 
       EXTERNAL ISMIN,ISMAX
@@ -30,5 +33,5 @@
 ! i.e based only on lat & lon values of each physical point.
 ! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48
-! For vizualisation : soon ...
+! For vizualisation : /u/tnalmd/bin/watercaps gcm_txt_output
        LOGICAL,SAVE :: improvedicelocation = .true.
 c
@@ -55,5 +58,9 @@
 
          alternate = 0
-         temptag = .true.
+         
+         if (ngridmx .ne. 1) then
+           watercaptag(:) = .false.
+           longwatercaptag(:) = .false.
+         endif
          
          write(*,*) "Ice dryness ?"
@@ -61,8 +68,10 @@
          call getin("icedryness",icedryness)
          write(*,*) " icedryness = ",icedryness
-
+         
         
 #ifdef MESOSCALE
+
       do ig=1,ngridmx
+
          !write(*,*) "all qsurf to zero. dirty."
          do iq=1,nqmx
@@ -81,39 +90,59 @@
                  watercaptag(ig)  = .false.
                  dryness(ig)      = 1.
-         endif  
-      enddo
+         endif 
+         
+      endif
 #else
-       ! print*,'ngridmx',ngridmx
 
       if (improvedicelocation) then
       
         if (ngridmx .eq. 738) then ! hopefully 32x24
+           
           print*,'water ice caps distribution for 32x24 resolution:'
-          watercaptag(1:9)    = .true. ! central cap - core
-          watercaptag(26:33)  = .true. ! central cap
+          longwatercaptag(1:9)    = .true. ! central cap - core
+          longwatercaptag(26:33)  = .true. ! central cap
             
         else if (ngridmx .eq. 3010) then ! hopefully 64x48
+
+! For latitudes:
+! [ 67.5   71.25  75.    78.75  82.5   86.25]
+! The water ice caps should be (according to TES temperatures):
+! [  8.63e-03   1.02e+00   5.99e+00   2.66e+01   5.65e+01]
+
           print*,'water ice caps distribution for 64x48 resolution:'
-          watercaptag(1:65)   = .true. ! central cap - core
-!          watercaptag(85)     = .true. ! central cap
-          watercaptag(83:85)  = .true. ! central cap - BIGGER
-          watercaptag(93:114) = .true. ! central cap
-!          watercaptag(254)    = .true. ! Korolev crater
-!          watercaptag(255)    = .true. ! Korolev crater --DECALE
-          watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG
-!--------------------------------------------------------------------          
-!          watercaptag(136)    = .true. ! outlier, lat = 78.75
-!          watercaptag(138)    = .true. ! outlier, lat = 78.75
-!          watercaptag(140)    = .true. ! outlier, lat = 78.75
-!          watercaptag(142)    = .true. ! outlier, lat = 78.75
-!          watercaptag(183)    = .true. ! outlier, lat = 78.75
-!          watercaptag(185)    = .true. ! outlier, lat = 78.75
-!          watercaptag(187)    = .true. ! outlier, lat = 78.75
-!          watercaptag(189)    = .true. ! outlier, lat = 78.75
-!          watercaptag(191)    = .true. ! outlier, lat = 78.75
-!          watercaptag(193)    = .true. ! outlier, lat = 78.75
-!--------------------------------------------------------------------          
-          watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75
-          watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75
+          longwatercaptag(1:65)   = .true. ! central cap - core
+          longwatercaptag(75:85)  = .true. ! central cap 
+          longwatercaptag(93:114) = .true. ! central cap
+!---------------------   OUTLIERS  ----------------------------          
+          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
+          longwatercaptag(194)    = .true. ! outlier, lat = 75
+          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
+!!-------------------------   OLD   ----------------------------          
+!!          longwatercaptag(83:85)  = .true. 
+!!          longwatercaptag(135:142)  = .true. 
+!!          longwatercaptag(181:193)  = .true. 
+!!          longwatercaptag(242:257)  = .true. 
+
             
         else if (ngridmx .ne. 1) then
@@ -129,7 +158,7 @@
         do ig=1,ngridmx
           dryness (ig) = icedryness
-          !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs
-          !!! OU remettre comme c'était avant Aymeric et pis c'est tout
-          if (watercaptag(ig)) then
+
+          if (longwatercaptag(ig)) then
+             watercaptag(ig) = .True.
              print*,'ice water cap', lati(ig)*180./pi,
      .              long(ig)*180./pi, ig
@@ -137,19 +166,5 @@
         enddo
 
-!!------------------------ test -----------------------------
-!!-----------------------------------------------------------    
-!!      else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA ....
-!!        
-!!         do ig=1,ngridmx
-!!          dryness (ig) = 1
-!!          if (albedodat(ig) .ge. 0.35) then
-!!             watercaptag(ig)    = .true.
-!!             print*,' albedo criteria ice water cap', lati(ig)*180./pi,
-!!     .              long(ig)*180./pi, ig
-!!          endif
-!!        enddo
-!!-----------------------------------------------------------    
-!!-----------------------------------------------------------    
-      
+
       else
 
@@ -213,23 +228,15 @@
      .         then
      
-             if (temptag) then
              
                if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
                 if (ngridmx.ne.1) watercaptag(ig)=.true.
-!                  print*,'ice water cap', lati(ig)*180./pi,
-!     .              long(ig)*180./pi, ig
-                  !dryness(ig) = 1.
+                  print*,'ice water cap', lati(ig)*180./pi,
+     .              long(ig)*180./pi, ig
+                  dryness(ig) = 1.
                   alternate = 1
                 else
                   alternate = 0
                 endif !end if alternate = 0
-               
-             else
-             
-!               if (ngridmx.ne.1) watercaptag(ig)=.true.
-!                  write(*,*) "outliers ", lati(ig)*180./pi,
-!     .              long(ig)*180./pi
-     
-             endif ! end if temptag
+
              
            endif
@@ -259,5 +266,5 @@
      .            long(ig)*180./pi, ig
               endif
-              !dryness(ig) = 1.
+             dryness(ig) = 1.
            endif
 
@@ -266,5 +273,5 @@
 c-------------------------------------------------------- 
 
-           if (lati(ig)*180./pi.gt.84) then
+           if (abs(lati(ig)*180./pi).gt.80) then
              print*,'ice water cap', lati(ig)*180./pi,
      .         long(ig)*180./pi, ig
Index: /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 633)
@@ -122,4 +122,5 @@
      &            ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
      &            (ccn0*4./3.*pi)),rdust(ig,l) )
+                rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
                 nuice(ig,l) = nuice_ref
               ENDDO
@@ -142,4 +143,5 @@
                 rice(ig,l) =
      &           CBRT( Mo/No * 0.75 / pi / rhocloud(ig,l))
+                rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
                 nuice(ig,l) = nuice_ref
               ENDDO
Index: /trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 632)
+++ /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 633)
@@ -1,8 +1,11 @@
-       SUBROUTINE watercloud(ngrid,nlay, ptimestep, 
-     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
+       SUBROUTINE watercloud(ngrid,nlay,ptimestep, 
+     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
      &                pq,pdq,pdqcloud,pdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud)
+! to use  'getin'
+      USE ioipsl_getincom
       IMPLICIT NONE
+
 
 c=======================================================================
@@ -13,8 +16,11 @@
 c    - An improved microphysical scheme (see improvedclouds.F)
 c
+c  There is a time loop specific to cloud formation 
+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 - Oct. 2011
+c  2004 - 2012
 c=======================================================================
 
@@ -35,27 +41,26 @@
 
       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)         ! tendance surf pressure
-      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
+      REAL pdpsrf(ngrid)         ! tendence surf pressure
       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)       ! tendance temperature des autres param.
+      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
 
       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
-      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,nlayermx) ! Dust geometric mean radius (m)
+      real pdq(ngrid,nlay,nq)    ! tendence 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)
 
 c   Outputs:
 c   -------
 
-      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 pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
+      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
+                                   ! a la chaleur latente
 
       REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
@@ -63,11 +68,34 @@
       REAL nuice(ngrid,nlay)   ! Estimated effective variance
                                !   of the size distribution
-      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
-      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
+      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
+      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
 
 c   local:
 c   ------
-
-      INTEGER ig,l
+      
+      ! for ice radius computation
+      REAL Mo,No
+      REAl ccntyp
+      
+      ! 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+physics)
+      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
+      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
+      
+      REAL CBRT
+      EXTERNAL CBRT
+
+
+      INTEGER iq,ig,l
       LOGICAL,SAVE :: firstcall=.true.
 
@@ -92,22 +120,279 @@
         write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
         write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
+                
+        write(*,*) "time subsampling for microphysic ?"
+#ifdef MESOSCALE
+        imicro = 2
+#else
+        imicro = 15
+#endif
+        call getin("imicro",imicro)
+        write(*,*)"imicro = ",imicro
+        
+        microtimestep = ptimestep/real(imicro)
+        write(*,*)"Physical timestep is",ptimestep 
+        write(*,*)"Microphysics timestep is",microtimestep 
 
         firstcall=.false.
       ENDIF ! of IF (firstcall)
       
-
-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,ptimestep,
-     &             pplev,pplay,pzlev,pzlay,pt,pdt,
-     &             pq,pdq,pdqcloud,pdtcloud,
-     &             nq,tau,rice,nuice,rsedcloud)
-      ENDIF
+c-----Initialization
+      subpdq(1:ngrid,1:nlay,1:nq) = 0
+      subpdt(1:ngrid,1:nlay)      = 0
+      
+      ! default value if no ice
+      rhocloud(1:ngrid,1:nlay) = rho_dust
+
+
+
+c------------------------------------------------------------------
+c Time subsampling for microphysics 
+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------ Tracers 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,
+     &             pplay,pt,subpdt,
+     &             pq,subpdq,subpdqcloud,subpdtcloud,
+     &             nq,tauscaling)
+
+        ELSE
+           CALL simpleclouds(ngrid,nlay,microtimestep,
+     &             pplay,pzlay,pt,subpdt,
+     &             pq,subpdq,subpdqcloud,subpdtcloud,
+     &             nq,tau)
+        ENDIF
+
+
+c-------------------------------------------------------------------
+c   3.  Updating tendencies after cloud scheme:
+c-----------------------------------------------
+
+        IF (microphys) THEN
+          DO l=1,nlay
+            DO ig=1,ngrid
+              subpdq(ig,l,igcm_dust_mass) =
+     &            subpdq(ig,l,igcm_dust_mass)
+     &          + subpdqcloud(ig,l,igcm_dust_mass)
+              subpdq(ig,l,igcm_dust_number) =
+     &            subpdq(ig,l,igcm_dust_number)
+     &          + subpdqcloud(ig,l,igcm_dust_number)
+              subpdq(ig,l,igcm_ccn_mass) =
+     &            subpdq(ig,l,igcm_ccn_mass)
+     &          + subpdqcloud(ig,l,igcm_ccn_mass)
+              subpdq(ig,l,igcm_ccn_number) =
+     &            subpdq(ig,l,igcm_ccn_number)
+     &          + subpdqcloud(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)
+     &        + subpdqcloud(ig,l,igcm_h2o_ice)
+            subpdq(ig,l,igcm_h2o_vap) =
+     &          subpdq(ig,l,igcm_h2o_vap)
+     &        + subpdqcloud(ig,l,igcm_h2o_vap)
+          ENDDO
+        ENDDO
+     
+ 
+      ENDDO ! of DO microstep=1,imicro
+      
+c-------------------------------------------------------------------
+c   6.  Compute final tendencies after time loop:
+c------------------------------------------------
+
+c------ Temperature tendency pdtcloud
+       DO l=1,nlay
+         DO ig=1,ngrid
+             pdtcloud(ig,l) =
+     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
+          ENDDO
+       ENDDO
+c------ Tracers tendencies pdqcloud
+       DO iq=1,nq
+        DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/real(imicro)
+     &       - pdq(ig,l,iq)
+         ENDDO
+        ENDDO
+       ENDDO
+
+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. 0.) THEN
+         pdqcloud(ig,l,igcm_ccn_number) =
+     &     - pq(ig,l,igcm_ccn_number)/ptimestep 
+     &     - pdq(ig,l,igcm_ccn_number)
+         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)
+         pdqcloud(ig,l,igcm_dust_mass) = 
+     &     -pdqcloud(ig,l,igcm_ccn_mass)
+          ENDIF
+         ENDDO
+        ENDDO
+       ENDIF
+
+        DO l=1,nlay
+         DO ig=1,ngrid
+          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
+     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice)) 
+     &       .le. 1.e-8) THEN
+           pdqcloud(ig,l,igcm_h2o_ice) = 
+     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
+           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
+          ENDIF 
+          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
+     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap)) 
+     &       .le. 1.e-8) THEN
+           pdqcloud(ig,l,igcm_h2o_vap) = 
+     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
+           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
+          ENDIF
+         ENDDO
+        ENDDO
+
+
+c------Update the ice and dust particle size "rice" for output or photochemistry
+c------Only rsedcloud is used for the water cycle
+
+      IF(scavenging) THEN 
+        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)
+
+          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
+        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)
+
+            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 
+
+
+!--------------------------------------------------------------
+!--------------------------------------------------------------
       
 
@@ -123,5 +408,6 @@
       end do
 
-      RETURN
+c=======================================================================
+
       END
  
