Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 728)
@@ -80,5 +80,5 @@
       REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
       ! BENJAMIN MODIFS
-      real CLFtot,CLF1,CLF2
+      real CLFtot
       real totcloudfrac(ngridmx)
       logical clearsky
@@ -145,16 +145,12 @@
 !==================================================================
 
-        if (iaero_co2.ne.0) then
+      if (iaero_co2.ne.0) then
            iaer=iaero_co2
 !       1. Initialization
-            aerosol(:,:,iaer)=0.0
+            aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
 !       2. Opacity calculation
             if (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
-               do ig=1, ngrid
-                  do l=1,nlayer
-                     aerosol(ig,l,iaer)=1.e-9
-                  end do
-                  !aerosol(ig,12,iaer)=4.0 ! single cloud layer option
-               end do
+               aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
+               !aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
             else
                DO ig=1, ngrid
@@ -179,20 +175,16 @@
                ENDDO
             end if ! if fixed or varying
-        end if ! if CO2 aerosols   
+      end if ! if CO2 aerosols   
 !==================================================================
 !     Water ice / liquid 
 !==================================================================
 
-        if (iaero_h2o.ne.0) then 
+      if (iaero_h2o.ne.0) then 
            iaer=iaero_h2o
 !       1. Initialization
-            aerosol(:,:,iaer)=0.0
+            aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
 !       2. Opacity calculation
             if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
-               do ig=1,ngrid ! to stop the rad tran bug
-                  do l=1,nlayer
-                     aerosol(ig,l,iaer) =1.e-9
-                  end do
-               end do
+               aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
 
                ! put cloud at cloudlvl
@@ -228,49 +220,51 @@
 
             else
+
                do ig=1, ngrid
                   do l=1,nlayer-1 ! to stop the rad tran bug
 
-                     if(CLFvarying)then
-                        CLF1    = max(cloudfrac(ig,l),0.01)
-                        CLFtot  = max(totcloudfrac(ig),0.01)
-                        CLF2    = CLF1/CLFtot 
-                        CLF2    = min(1.,CLF2)
-                     else
-                        CLF1=1.0
-                        CLF2=CLFfixval
-                     endif
-                     
-                     aerosol0 =                                   &
+                     aerosol(ig,l,iaer) =                                    & !modification by BC
                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
                           ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
                           !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
-                          !( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the 
-                          !  CLF1				    !   clear=false/pq=0 case
+                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the 
+                                                                     !   clear=false/pq=0 case
                           ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   &  
-                            CLF1
-
-                     aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0)) 
-                     ! why no division in exponential?
-                     ! because its already done in aerosol0
-                     
-                     aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
-                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)                     
-                     
+                          ( pplev(ig,l) - pplev(ig,l+1) ) / g 
+
                   enddo
                enddo
-            end if
-        end if ! End if h2o aerosol
+
+               if(CLFvarying)then
+                  call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
+                  do ig=1, ngrid
+                     do l=1,nlayer-1 ! to stop the rad tran bug
+                        CLFtot  = max(totcloudfrac(ig),0.01)
+                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
+                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
+                     enddo
+                  enddo
+               else
+                  do ig=1, ngrid
+                     do l=1,nlayer-1 ! to stop the rad tran bug
+                        CLFtot  = CLFfixval
+                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
+                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
+                     enddo
+                  enddo
+              end if!(CLFvarying)
+            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
+	      
+      end if ! End if h2o aerosol
 
 !==================================================================
 !             Dust 
 !==================================================================
-        if (iaero_dust.ne.0) then
+      if (iaero_dust.ne.0) then
           iaer=iaero_dust
-
 !         1. Initialization 
-            aerosol(:,:,iaer)=0.0
+          aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
           
-            topdust=30.0 ! km  (used to be 10.0 km) LK
+          topdust=30.0 ! km  (used to be 10.0 km) LK
 
 !       2. Opacity calculation
@@ -314,5 +308,5 @@
               ENDDO
             ENDDO
-        end if ! If dust aerosol   
+      end if ! If dust aerosol   
 
 !==================================================================
@@ -321,8 +315,8 @@
 ! added by LK
       if (iaero_h2so4.ne.0) then
-          iaer=iaero_h2so4
+         iaer=iaero_h2so4
 
 !       1. Initialization
-            aerosol(:,:,iaer)=0.0
+         aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
 
 
@@ -330,28 +324,25 @@
 
 !           expfactor=0.
-           DO l=1,nlayer-1
-             DO ig=1,ngrid
-!             Typical mixing ratio profile
-
-                 zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
-                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
+         DO l=1,nlayer-1
+            DO ig=1,ngrid
+!              Typical mixing ratio profile
+
+               zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
+               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
 
 !             Vertical scaling function
-              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
-               *expfactor
-
-
-             ENDDO
-           ENDDO
-            tauh2so4tmp(1:ngrid)=0.
-              DO l=1,nlayer
-                DO ig=1,ngrid
-                   tauh2so4tmp(ig) = tauh2so4tmp(ig) &
-                          +  aerosol(ig,l,iaer)
-                ENDDO
-              ENDDO
-            DO l=1,nlayer-1
-               DO ig=1,ngrid
-                  aerosol(ig,l,iaer) = max(1E-20, &
+               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
+
+            ENDDO
+         ENDDO
+         tauh2so4tmp(1:ngrid)=0.
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
+            ENDDO
+         ENDDO
+         DO l=1,nlayer-1
+            DO ig=1,ngrid
+               aerosol(ig,l,iaer) = max(1E-20, &
                           1 &
                        *  pplev(ig,1) / preff &
@@ -359,6 +350,6 @@
                        /  tauh2so4tmp(ig))
 
-              ENDDO
             ENDDO
+         ENDDO
 
 ! 1/700. is assuming a "sulfurtau" of 1
Index: trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 728)
@@ -5,8 +5,8 @@
 !  aerosol indexes: these are initialized in iniaerosol.F and should be 0 if the
 !                 corresponding aerosol was not activated in callphys.def
-      integer :: iaero_co2 
-      integer :: iaero_h2o
-      integer :: iaero_dust
-      integer :: iaero_h2so4
+      integer, save :: iaero_co2 
+      integer, save :: iaero_h2o
+      integer, save :: iaero_dust
+      integer, save :: iaero_h2so4
 
 !==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 728)
@@ -14,4 +14,5 @@
       use ioipsl_getincom 
       use gases_h
+      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
       use aerosol_mod
 
@@ -74,5 +75,5 @@
 
       REAL reffrad(ngrid,nlayer,naerkind)
-      REAL nueffrad(ngrid,nlayer,naerkind)
+      REAL, SAVE :: nueffrad(ngridmx,nlayermx,naerkind)
 
 !     OUTPUT
@@ -171,15 +172,6 @@
       real pqtest(ngridmx,nlayer,nq)
 
-!     are particle radii fixed?
-      logical, save ::  radfixed
       real maxrad, minrad
-
-!     Local water cloud optical properties
-      real, parameter ::  rad_froid=35.0e-6 
-      real, parameter ::  rad_chaud=10.0e-6 
-      real, parameter ::  coef_chaud=0.13
-      real, parameter ::  coef_froid=0.09
-      real zfice
-      
+            
       real CBRT
       external CBRT
@@ -201,10 +193,6 @@
       qsiaer(:,:,:)=0.0
       giaer(:,:,:) =0.0
-      radfixed=.false.
 
       if(firstcall) then
-         if(kastprof)then
-            radfixed=.true.
-         endif  
 
          call system('rm -f surf_vals_long.out')
@@ -212,46 +200,14 @@
 !--------------------------------------------------
 !     Effective radius and variance of the aerosols
-
-
-!     these values will change once the microphysics gets to work
-!     UNLESS tracer=.false., in which case we should be working with
-!     a fixed aerosol layer, and be able to define reffrad in a 
-!     .def file. To be improved!
-!     Generalization of aerosols added by LK
-            do iaer = 1,naerkind
-              if (iaer.eq.iaero_co2) then ! active CO2 aerosols  
-                do l=1,nlayer
-                   do ig=1,ngrid
-                     reffrad(ig,l,iaer)  = 1.e-4
-                     nueffrad(ig,l,iaer) = 0.1 
-                   enddo
-                enddo
-              endif
-              if (iaer.eq.iaero_h2o) then ! active H2O aerosols
-                do l=1,nlayer
-                   do ig=1,ngrid
-                     reffrad(ig,l,iaer)  = 1.e-5
-                     nueffrad(ig,l,iaer) = 0.1 
-                   enddo
-                enddo
-              endif
-              if(iaer.eq.iaero_dust)then ! active dust
-                do l=1,nlayer
-                   do ig=1,ngrid
-                     reffrad(ig,l,iaer)  = 1.e-5
-                     nueffrad(ig,l,iaer) = 0.1 
-                   enddo
-                enddo
-              endif
-              if(iaer.eq.iaero_h2so4)then ! Active h2so4 aerosols
-                do l=1,nlayer
-                   do ig=1,ngrid
-                     reffrad(ig,l,iaer)  = 1.e-6
-                     nueffrad(ig,l,iaer) = 0.1
-                   enddo
-                enddo
-              endif
-            enddo
-
+         if(naerkind.gt.4)then
+            print*,'Code not general enough to deal with naerkind > 4 yet.'
+            call abort
+         endif
+         call su_aer_radii(reffrad,nueffrad)
+	 
+	 
+
+!--------------------------------------------------
+!     set up correlated k
          print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
          call getin("corrkdir",corrkdir)
@@ -305,112 +261,27 @@
 !     Initialization on every call    
 
-      do l=1,nlayer
-         do ig=1,ngrid
-            do iaer=1,naerkind
-               nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
-            enddo
-         enddo
-      enddo
-
-
+!--------------------------------------------------
+!     Effective radius and variance of the aerosols
       do iaer=1,naerkind
 
-       if(iaer.eq.iaero_co2)then
-         if(radfixed)then
-           do l=1,nlayer
-            do ig=1,ngrid
-               reffrad(ig,l,iaer)    = 5.e-5 ! CO2 ice
-            enddo
-           enddo
-           print*,'CO2 ice particle size =',reffrad(1,1,iaer)/1.e-6,'um'
-         else 
-           maxrad=0.0
-           !averad=0.0
-           minrad=1.0
-           do l=1,nlayer
-
-             !masse = (pplev(ig,l) - pplev(ig,l+1))/g
-
-            do ig=1,ngrid
-               if(tracer.and.igcm_co2_ice.gt.0)then
-
-                  if(igcm_co2_ice.lt.1)then
-                     print*,'Tracers but no CO2 ice still'
-                     print*,'seems to be a problem...'
-                     print*,'Aborting in callcorrk.'
-                     stop
-                  endif
-                  reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
-                       (4*Nmix_co2*pi*rho_co2) )
-               endif
-               reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
-               reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),500.e-6)
-               !averad = averad + reffrad(ig,l,iaer)*area(ig)
-               maxrad = max(reffrad(ig,l,iaer),maxrad)
-               minrad = min(reffrad(ig,l,iaer),minrad)
-            enddo
-           enddo
-           if(igcm_co2_ice.gt.0)then
-             print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
-             print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
-           endif
-         endif
-       endif
-
-
-       if(iaer.eq.iaero_h2o)then
-          if(radfixed) then
-            do l=1,nlayer
-               do ig=1,ngrid
-                  !reffrad(ig,l,iaer) = 2.e-5 ! H2O ice
-                  reffrad(ig,l,iaer) = rad_chaud ! H2O ice
-
-                  zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
-                  zfice = MIN(MAX(zfice,0.0),1.0)
-                  reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
-                  nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
-               enddo
-            enddo
-            print*,'H2O ice particle size=',reffrad(1,1,iaer)/1.e-6,'um'
-          else
-             if(water)then
-               maxrad=0.0
-               minrad=1.0
-               do l=1,nlayer
-                do ig=1,ngrid
-                  reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
-                       (4*Nmix_h2o*pi*rho_ice) )
-                  reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
-                  reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),100.e-6)
-
-                  maxrad = max(reffrad(ig,l,iaer),maxrad)
-                  minrad = min(reffrad(ig,l,iaer),minrad)
-                enddo
-               enddo
-               print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
-               print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
-             endif
-         endif
-       endif
-
+         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
+	    call co2_reffrad(pq,reffrad)
+            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
+            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
+	 end if
+         if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
+	    call h2o_reffrad(pq,pt,reffrad,nueffrad)
+            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
+            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
+         endif
          if(iaer.eq.iaero_dust)then
-            do l=1,nlayer
-               do ig=1,ngrid
-                  reffrad(ig,l,iaer) = 2.e-6 ! dust
-               enddo
-            enddo
+	    call dust_reffrad(reffrad)
             print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
-
          if(iaer.eq.iaero_h2so4)then
-            do l=1,nlayer
-               do ig=1,ngrid
-                  reffrad(ig,l,iaer) = 1.e-6 ! h2so4
-               enddo
-            enddo
+	    call h2so4_reffrad(reffrad)
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
-      enddo ! do iaer=1,naerkind
-         ! added by LK
+      end do !iaer=1,naerkind
 
 
@@ -641,7 +512,7 @@
       if(.not.kastprof)then
       ! IMPORTANT: Now convert from kg/kg to mol/mol
-      do k=1,L_LEVELS
-         qvar(k) = qvar(k)/epsi
-      end do
+         do k=1,L_LEVELS
+            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
+         end do
       end if
 
@@ -764,6 +635,8 @@
             print*,'Maximum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
-            !print*,'WARNING, OVERRIDING FOR TEST'
-            call abort
+	    print*,'level,grid,T',k,ig,tlevrad(k)
+            print*,'WARNING, OVERRIDING FOR TEST'
+            !call abort
+	    tlevrad(k)=tgasmax
          endif
       enddo
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 728)
@@ -9,39 +9,21 @@
      &   , enertest                                                     &
      &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
-     &   , Nmix_co2, Nmix_h2o                                           &
-     &   , dusttau                                                      &
-     &   , meanOLR                                                      &
-     &   , specOLR                                                      &
-     &   , kastprof                                                     &
-     &   , noradsurf                                                    &
-     &   , Tstrat                                                       &
-     &   , newtonian                                                    &
-     &   , check_cpp_match                                              &
-     &   , force_cpp                                                    &
-     &   , tau_relax                                                    &
-     &   , testradtimes                                                 &
+     &   , Nmix_co2, radfixed, dusttau                                  &
+     &   , meanOLR, specOLR                                             &
+     &   , kastprof, noradsurf, Tstrat                                  &
+     &   , newtonian, tau_relax, testradtimes                           &
+     &   , check_cpp_match, force_cpp                                   &
      &   , rayleigh                                                     &
-     &   , stelbbody                                                    &
-     &   , stelTbb                                                      &
+     &   , stelbbody, stelTbb                                           &
      &   , tplanet                                                      &
-     &   , startype                                                     &
-     &   , Fat1AU                                                       &
+     &   , startype, Fat1AU                                             &
      &   , nearco2cond                                                  &
-     &   , tracer                                                       &
-     &   , varactive                                                    &
-     &   , varfixed                                                     &
-     &   , satval                                                       &
+     &   , tracer, mass_redistrib, varactive, varfixed, satval          &
      &   , sedimentation,water,watercond,waterrain                      &
-     &   , rainthreshold                                                &
      &   , aeroco2,aeroh2o,aeroh2so4                                    &
      &   , aerofixco2,aerofixh2o                                        &
-     &   , hydrology                                                    &
-     &   , sourceevol                                                   &
-     &   , icetstep                                                     &
-     &   , albedosnow                                                   &
-     &   , maxicethick                                                  &
-     &   , Tsaldiff                                                     &
-     &   , CLFfixval                                                    &
-     &   , CLFvarying                                                   &
+     &   , hydrology, sourceevol, icetstep, albedosnow                  &
+     &   , maxicethick, Tsaldiff                                        &
+     &   , CLFfixval, CLFvarying                                        &
      &   , n2mixratio                                                   &
      &   , co2supsat                                                    &
@@ -69,6 +51,8 @@
       logical nearco2cond
       logical tracer
+      logical mass_redistrib
       logical varactive
       logical varfixed
+      logical radfixed
       logical sedimentation
       logical water,watercond,waterrain
@@ -87,5 +71,4 @@
       real topdustref
       real Nmix_co2
-      real Nmix_h2o
       real dusttau
       real Fat1AU
@@ -94,5 +77,4 @@
       real tplanet
       real satval
-      real rainthreshold
       real CLFfixval
       real n2mixratio
Index: trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 728)
@@ -8,4 +8,6 @@
       use radinc_h, only : naerkind
       use gases_h
+      use radii_mod, only : co2_reffrad
+      use aerosol_mod, only : iaero_co2
 
       implicit none
@@ -146,6 +148,4 @@
       character(len=3) :: gasname
 
-      real reffradmin, reffradmax
-
       real ppco2
 
@@ -157,12 +157,4 @@
       call zerophys(ngridmx*nlayermx*nqmx,zq)
       call zerophys(ngridmx*nlayermx,zt)
-
-      !reffradmin=1.e-7
-      !reffradmax=5.e-4
-      !reffradmin=0.5e-7
-      !reffradmax=0.1e-3 ! FF data
-      reffradmin=0.1e-5
-      reffradmax=0.1e-3 ! JB data
-!     improve this in the future...
 
 !     Initialisation
@@ -319,24 +311,11 @@
 !     Gravitational sedimentation
             
-!     sedimentation computed from radius computed from q
-!     assuming that the ice is splitted in Nmix particle
+!     sedimentation computed from radius computed from q in module radii_mod
+	 call co2_reffrad(zq,reffrad)
+	 
          do  ilay=1,nlayer
             do ig=1, ngrid
 
-               reff = CBRT( 3*zq(ig,ilay,i_co2ice)/( 4*Nmix_co2*pi*rho_co2 ))
-
-               ! there should be a more elegant way of doing this...
-               if(reff.lt.1.e-16) reff=1.e-16
-
-               ! update reffrad for radiative transfer
-               if(reff.lt.reffradmin)then
-                  reffrad(ig,ilay,1)=reffradmin
-                  !print*,'reff below optical limit'
-               elseif(reff.gt.reffradmax)then
-                  reffrad(ig,ilay,1)=reffradmax
-                  !print*,'reff above optical limit'
-               else
-                  reffrad(ig,ilay,1)=reff
-               endif
+               reff = reffrad(ig,ilay,iaero_co2)
 
                call stokes                      &
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 728)
@@ -132,4 +132,10 @@
          write(*,*) " tracer = ",tracer
 
+         write(*,*) "Run with or without atm mass update ",
+     &      " due to tracer evaporation/condensation?"
+         mass_redistrib=.false. ! default value
+         call getin("mass_redistrib",mass_redistrib)
+         write(*,*) " mass_redistrib = ",mass_redistrib
+
          write(*,*) "Diurnal cycle ?"
          write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
@@ -377,13 +383,17 @@
          write(*,*)" CLFfixval = ",CLFfixval
 
-         write(*,*)"Number mixing ratio of CO2 ice particles:"
-         Nmix_co2=100000. ! default value
+         write(*,*)"fixed radii for Cloud particles?"
+         radfixed=.false. ! default value
+         call getin("radfixed",radfixed)
+         write(*,*)" radfixed = ",radfixed
+
+         if(kastprof)then
+            radfixed=.true.
+         endif  
+
+	 write(*,*)"Number mixing ratio of CO2 ice particles:"
+         Nmix_co2=1.e6 ! default value
          call getin("Nmix_co2",Nmix_co2)
          write(*,*)" Nmix_co2 = ",Nmix_co2
-
-         write(*,*)"Number mixing ratio of H2O ice particles:"
-         Nmix_h2o=10000000. ! default value
-         call getin("Nmix_h2o",Nmix_h2o)
-         write(*,*)" Nmix_h2o = ",Nmix_h2o
 
 !         write(*,*)"Number of radiatively active aerosols:"
@@ -499,9 +509,4 @@
          call getin("waterrain",waterrain)
          write(*,*) " waterrain = ",waterrain
-
-         write(*,*) "Precipitation threshold ?"
-         rainthreshold=0.011 ! default value (Emmanuel 1997)
-         call getin("rainthreshold",rainthreshold)
-         write(*,*) " rainthreshold = ",rainthreshold
 
          write(*,*) "Include surface hydrology ?"
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 728)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 728)
@@ -0,0 +1,148 @@
+      subroutine largescale(ptimestep, pplev, pplay, pt, pq,   &
+                        pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
+
+      use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
+          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
+      IMPLICIT none
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Calculates large-scale (stratiform) H2O condensation.
+!     
+!     Authors
+!     -------
+!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
+!     Original author Z. X. Li (1993)
+!     
+!==================================================================
+
+#include "dimensions.h" 
+#include "dimphys.h"
+#include "comcstfi.h"
+
+#include "fisice.h"
+#include "callkeys.h"
+#include "tracer.h"
+
+
+
+!     Arguments
+      REAL ptimestep                 ! intervalle du temps (s)
+      REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche
+      REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche
+      REAL pt(ngridmx,nlayermx)      ! temperature (K)
+      real pq(ngridmx,nlayermx,nqmx) ! tracer mixing ratio (kg/kg)
+      REAL pdt(ngridmx,nlayermx)     ! physical temperature tenedency (K/s)
+      REAL pdq(ngridmx,nlayermx,nqmx)! physical tracer tenedency (K/s)
+      REAL pdtlsc(ngridmx,nlayermx)  ! incrementation de la temperature (K)
+      REAL pdqvaplsc(ngridmx,nlayermx) ! incrementation de la vapeur d'eau
+      REAL pdqliqlsc(ngridmx,nlayermx) ! incrementation de l'eau liquide
+      REAL rneb(ngridmx,nlayermx)    ! fraction nuageuse
+
+
+!     Options du programme
+      REAL ratqs   ! determine largeur de la distribution de vapeur
+      PARAMETER (ratqs=0.2)
+
+!     Variables locales
+      REAL CBRT
+      EXTERNAL CBRT
+      INTEGER i, k , nn
+      INTEGER,PARAMETER :: niter=4
+      REAL zt(ngridmx), zq(ngridmx)
+      REAL zcond(ngridmx),zcond_iter
+      REAL zdelq(ngridmx)
+      REAL zqs(ngridmx), zdqs(ngridmx)
+      REAL psat_tmp
+      
+! evaporation calculations
+      REAL dqevap(ngridmx,nlayermx),dtevap(ngridmx,nlayermx)     
+      REAL qevap(ngridmx,nlayermx,nqmx)
+      REAL tevap(ngridmx,nlayermx)
+
+      REAL zcor(ngridmx), zdelta(ngridmx), zcvm5(ngridmx)
+      REAL zx_q(ngridmx)
+      REAL Nmix_local,zfice
+
+!     GCM -----> subroutine variables, initialisation of outputs
+
+      pdtlsc(1:ngridmx,1:nlayermx)  = 0.0
+      pdqvaplsc(1:ngridmx,1:nlayermx)  = 0.0
+      pdqliqlsc(1:ngridmx,1:nlayermx) = 0.0
+      rneb(1:ngridmx,1:nlayermx) = 0.0
+
+
+      ! Evaporate cloud water/ice
+      call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
+      ! note: we use qevap but not tevap in largescale/moistadj
+            ! otherwise is a big mess
+
+
+!  Boucle verticale (du haut vers le bas)
+   DO k = nlayermx, 1, -1
+
+      zt(1:ngridmx)=pt(1:ngridmx,k)+(pdt(1:ngridmx,k)+dtevap(1:ngridmx,k))*ptimestep
+      zq(1:ngridmx)=qevap(1:ngridmx,k,igcm_h2o_vap) !liquid water is included in qevap
+
+!     Calculer la vapeur d'eau saturante et 
+!     determiner la condensation partielle
+      DO i = 1, ngridmx
+
+         if(zt(i).le.15.) zt(i)=15.   ! check too low temperatures
+!         call watersat(zt(i),pplay(i,k),zqs(i))
+         call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+ 
+         zdelq(i) = ratqs * zq(i)
+	 if(zq(i)+zdelq(i).lt.0.999) then
+	    rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
+	    zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
+            if (rneb(i,k) .LE. 0.0) zx_q(i) = 0.0
+            if (rneb(i,k) .GE. 1.0) zx_q(i) = zq(i)
+            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
+	 else
+	    if(zq(i).gt.zqs(i)) then
+	       rneb(i,k)=1.
+	       zx_q(i)=zq(i)
+	    else
+	       rneb(i,k)=0.
+	       zx_q(i)=zqs(i)  !no condensation needed
+	    Endif
+	 Endif
+
+!        iterative process to stabilize the scheme when large water amounts JL12
+         zcond(i) = 0.0
+         Do nn=1,niter  
+!            call watersat_grad(zt(i),zqs(i),zdqs(i))
+	    call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+            zcond_iter = MAX(0.0,(zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)))	   
+               !zcond always postive! cannot evaporate clouds!
+               !this is why we must reevaporate before largescale
+            zx_q(i) = zx_q(i) - zcond_iter
+	    zcond(i) = zcond(i) + zcond_iter
+	    zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
+            call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+	 End do ! niter
+
+
+         zcond(i) = zcond(i)/ptimestep ! added by RDW
+      ENDDO
+
+!     Tendances de t et q
+         pdqvaplsc(1:ngridmx,k)  = dqevap(1:ngridmx,k) - zcond(1:ngridmx)
+         pdqliqlsc(1:ngridmx,k) = - pdqvaplsc(1:ngridmx,k)
+         pdtlsc(1:ngridmx,k)  = pdqliqlsc(1:ngridmx,k)*RLVTT/RCPD
+
+   Enddo ! k= nlayermx, 1, -1
+   
+      !print*,'qsat=',zqs
+      !print*,'q=',q
+      !print*,'dq=',pdqvaplsc*ptimestep
+      !print*,'dT in LS=',pdtlsc*ptimestep
+
+      !print*,'rice=',rice
+      !print*,'rneb=',rneb
+
+      return
+      end
Index: trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 728)
+++ trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 728)
@@ -0,0 +1,404 @@
+      SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
+                       rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
+		       pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
+		       pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
+                                                   
+       USE watercommon_h, Only: Tsat_water,RLVTT
+
+       IMPLICIT NONE
+!=======================================================================
+!   subject:
+!   --------
+!     Mass and momentum fluxes through sigma levels as the surface pressure is modified are also taken into account
+!
+!   author:   Jeremy Leconte 2012 (from F.Forget 1998)
+!   ------
+!
+!   input:
+!   ----- 
+!    ngrid                 nombre de points de verticales
+!                          (toutes les boucles de la physique sont au
+!                          moins vectorisees sur ngrid)
+!    nlayer                nombre de couches
+!    pplay(ngrid,nlayer)   Pressure levels 
+!    pplev(ngrid,nlayer+1) Pressure levels 
+!    nq                    Number of tracers
+!
+!    pt(ngrid,nlayer)      temperature (en K)
+!    pq(ngrid,nlayer,nq)   tracer specific concentration (kg/kg of air)
+!    pu,pv (ngrid,nlayer)  wind velocity (m/s)
+!
+!                    
+!    pdX                   physical tendency of X before mass redistribution
+!
+!    pdmassmr                air Mass added to the atmosphere in each layer (kg/m2/s)
+!
+!   output:
+!   -------
+!
+!    pdXmr(ngrid)           physical tendency of X after mass redistribution
+!
+!    
+!
+!=======================================================================
+!
+!    0.  Declarations :
+!    ------------------
+!
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "surfdat.h"
+#include "comgeomfi.h"
+#include "comvert.h"
+#include "paramet.h"
+#include "callkeys.h"
+#include "tracer.h"
+
+!-----------------------------------------------------------------------
+!    Arguments :
+!    ---------
+      INTEGER ngrid, nlayer, nq   
+      REAL ptimestep
+      REAL pcapcal(ngridmx)
+      INTEGER rnat(ngridmx)      
+      REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
+      REAL pt(ngrid,nlayer),pdt(ngrid,nlayer)
+      REAL ptsrf(ngrid),pdtsrf(ngrid)
+      REAL pdtmr(ngrid,nlayer)
+      REAL pu(ngrid,nlayer) ,  pv(ngrid,nlayer)
+      REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
+      REAL pdmassmr(ngrid,nlayer)
+      REAL pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer)
+      REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
+      REAL pqs(ngridmx,nq)
+      REAL pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq)
+      REAL pdpsrfmr(ngrid),pdtsrfmr(ngrid)
+!
+!    Local variables :
+!    -----------------
+
+!    Boiling/sublimation
+      REAL Tsat(ngridmx),zmassboil(ngridmx)
+
+!    vertical reorganization of sigma levels
+      REAL zzu(nlayermx),zzv(nlayermx)
+      REAL zzq(nlayermx,nq),zzt(nlayermx)
+!    Dummy variables      
+      INTEGER n,l,ig,iq
+      REAL zdtsig(ngridmx,nlayermx)
+      REAL zmass(ngridmx,nlayermx),zzmass(nlayermx),w(nlayermx+1)
+      REAL zdmass_sum(ngridmx,nlayermx+1)
+      REAL zmflux(nlayermx+1)
+      REAL zq1(nlayermx)
+      REAL ztsrf(ngridmx)
+      REAL ztm(nlayermx+1) 
+      REAL zum(nlayermx+1) , zvm(nlayermx+1)
+      REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
+
+!   local saved variables
+      LOGICAL, SAVE :: firstcall=.true.
+
+!----------------------------------------------------------------------
+
+!   Initialisation
+!   --------------
+!
+      IF (firstcall) THEN
+         firstcall=.false.	 
+      ENDIF
+!
+!======================================================================
+!    Calcul of h2o condensation  
+!    ============================================================
+!  
+!    Used variable :
+!       pdmassmr      : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
+!       zdmass_sum(ngrid,l) : total air mass added to the atm above layer l per unit time (kg/m2/s)
+!
+!
+!     Surface tracer Tendencies set to 0
+!     -------------------------------------
+      pdqsmr(1:ngridmx,1:nqmx)=0.
+
+      ztsrf(1:ngridmx) = ptsrf(1:ngridmx) + pdtsrf(1:ngridmx)*ptimestep
+
+
+      DO ig=1,ngrid
+         zdmass_sum(ig,nlayermx+1)=0.
+         DO l = nlayermx, 1, -1
+           zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/g
+	   zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
+         END DO
+      END DO
+
+
+      if (water) then
+         do ig=1,ngridmx
+	    call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)*g*ptimestep,Tsat(ig))
+	 enddo
+               call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
+	 
+         do ig=1,ngridmx
+	    if (ztsrf(ig).gt.Tsat(ig)) then
+	       zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep
+	       if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_h2o_vap)).and.(rnat(ig).eq.1)) then
+	          zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
+	       endif
+	       zmassboil(ig)=zmassboil(ig)*0.3 !momentary, should be 1. JL12
+               pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
+	       pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)
+	       ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep
+	    else
+	       zmassboil(ig)=0.
+	       pdtsrfmr(ig)=0.
+	    endif
+	 enddo
+      endif
+
+!     *************************
+!           UPDATE SURFACE
+!     *************************
+!    Changing pressure at the surface:
+!    """"""""""""""""""""""""""""""""""""
+         
+      pdpsrfmr(1:ngridmx) = (zdmass_sum(1:ngridmx,1)+zmassboil(1:ngridmx))*g
+
+      do ig = 1, ngridmx
+        IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN
+         PRINT*,'STOP in condens'
+         PRINT*,'condensing more than total mass'
+         PRINT*,'Grid point ',ig
+         PRINT*,'Ps = ',pplev(ig,1)
+         PRINT*,'d Ps = ',pdpsrfmr(ig)*ptimestep
+         STOP
+        ENDIF
+      enddo ! of DO ig=1,ngrid
+
+
+! ***************************************************************
+!  Correction to account for redistribution between sigma or hybrid 
+!  layers when changing surface pressure
+!  zzx quantities have dimension (nlayermx) to speed up calculation
+! *************************************************************
+
+      DO ig=1,ngrid
+         zzt(1:nlayermx)  = pt(ig,1:nlayermx) + pdt(ig,1:nlayermx) * ptimestep
+         zzu(1:nlayermx)  = pu(ig,1:nlayermx) + pdu(ig,1:nlayermx) * ptimestep
+         zzv(1:nlayermx)  = pv(ig,1:nlayermx) + pdv(ig,1:nlayermx) * ptimestep
+         zzq(1:nlayermx,1:nqmx)=pq(ig,1:nlayermx,1:nqmx)+pdq(ig,1:nlayermx,1:nqmx)*ptimestep ! must add the water that has fallen???
+
+!  Mass fluxes of air through the sigma levels (kg.m-2.s-1)  (>0 when up)
+!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+         zmflux(1) = zmassboil(ig)
+         DO l=1,nlayer
+            zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
+! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
+            if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
+         END DO
+
+! Mass of each layer
+! ------------------ 
+         zzmass(1:nlayermx)=zmass(ig,1:nlayermx)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
+
+
+!  Corresponding fluxes for T,U,V,Q
+!  """"""""""""""""""""""""""""""""
+
+!        averaging operator for TRANSPORT  
+!        """"""""""""""""""""""""""""""""
+!        Value transfert at the surface interface when condensation
+!        sublimation:
+         ztm(1) = ztsrf(ig)
+         zum(1) = 0.  
+         zvm(1) = 0.  
+         zqm(1,1:nqmx)=0. ! most tracer do not condense !
+         if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
+	 
+!        Van Leer scheme:
+         w(1:nlayermx+1)=-zmflux(1:nlayermx+1)*ptimestep
+         call vl1d(zzt,2.,zzmass,w,ztm) 
+         call vl1d(zzu ,2.,zzmass,w,zum) 
+         call vl1d(zzv ,2.,zzmass,w,zvm) 
+         do iq=1,nqmx
+           zq1(1:nlayermx)=zzq(1:nlayermx,iq)
+           zqm1(1)=zqm(1,iq)
+!		print*,iq
+!		print*,zq1
+           call vl1d(zq1,2.,zzmass,w,zqm1)
+           do l=2,nlayer
+              zzq(l,iq)=zq1(l)
+              zqm(l,iq)=zqm1(l)
+           enddo
+         enddo
+
+!        Surface condensation affects low winds
+         if (zmflux(1).lt.0) then 
+            zum(1)= zzu(1) *  (w(1)/zzmass(1))
+            zvm(1)= zzv(1) *  (w(1)/zzmass(1))
+            if (w(1).gt.zzmass(1)) then ! ensure numerical stability
+               zum(1)= (zzu(1)-zum(2))*zzmass(1)/w(1) + zum(2)
+               zvm(1)= (zzv(1)-zvm(2))*zzmass(1)/w(1) + zvm(2)
+            end if
+         end if
+                    
+         ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... 
+         zum(nlayer+1)= zzu(nlayer)  ! should not be used, but...
+         zvm(nlayer+1)= zzv(nlayer)  ! should not be used, but...
+         zqm(nlayer+1,1:nqmx)= zzq(nlayer,1:nqmx)
+ 
+!        Tendencies on T, U, V, Q 
+!        """"""""""""""""""""""""
+         DO l=1,nlayer
+ 
+!           Tendencies on T
+            pdtmr(ig,l) = (1/zzmass(l)) *   &
+		(zmflux(l)*(ztm(l) - zzt(l))-zmflux(l+1)*(ztm(l+1)-zzt(l)))
+		  !JL12 the last term in Newcondens has been set to zero because we are only dealing with redistribution here
+
+!           Tendencies on U
+            pdumr(ig,l)   = (1/zzmass(l)) *( zmflux(l)*(zum(l) - zzu(l)) - zmflux(l+1)*(zum(l+1) - zzu(l)) )
+
+!           Tendencies on V
+            pdvmr(ig,l)   = (1/zzmass(l)) *( zmflux(l)*(zvm(l) - zzv(l)) - zmflux(l+1)*(zvm(l+1) - zzv(l)) )
+
+         END DO
+
+!        Tendencies on Q
+         do iq=1,nqmx
+            DO l=1,nlayer
+               pdqmr(ig,l,iq)= (1/zzmass(l)) *   &
+		   (zmflux(l)*(zqm(l,iq)-zzq(l,iq))- zmflux(l+1)*(zqm(l+1,iq)-zzq(l,iq)) - pdmassmr(ig,l)*zzq(l,iq))
+            END DO
+         enddo
+
+      END DO  ! loop on ig 
+
+      return
+      end
+
+
+
+! *****************************************************************
+      SUBROUTINE vl1d(q,pente_max,zzmass,w,qm)
+!
+!    
+!     Operateur de moyenne inter-couche pour calcul de transport type
+!     Van-Leer " pseudo amont " dans la verticale
+!    q,w sont des arguments d'entree  pour le s-pg ....
+!    masse : masse de la couche Dp/g
+!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
+!    pente_max = 2 conseillee
+!
+!
+!   --------------------------------------------------------------------
+      IMPLICIT NONE
+
+#include "dimensions.h"
+
+!   Arguments:
+!   ----------
+      real zzmass(llm),pente_max
+      REAL q(llm),qm(llm+1)
+      REAL w(llm+1)
+!
+!      Local 
+!   ---------
+!
+      INTEGER l
+!
+      real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
+      real sigw, Mtot, MQtot
+      integer m 
+!     integer ismax,ismin 
+
+
+!    On oriente tout dans le sens de la pression 
+!     W > 0 WHEN DOWN !!!!!!!!!!!!!
+
+      do l=2,llm
+            dzqw(l)=q(l-1)-q(l)
+            adzqw(l)=abs(dzqw(l))
+      enddo
+
+      do l=2,llm-1
+            if(dzqw(l)*dzqw(l+1).gt.0.) then
+                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
+            else
+                dzq(l)=0.
+            endif
+            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
+            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
+      enddo
+
+         dzq(1)=0.
+         dzq(llm)=0.
+
+       do l = 1,llm-1
+
+!         Regular scheme (transfered mass < layer mass)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          if(w(l+1).gt.0. .and. w(l+1).le.zzmass(l+1)) then
+             sigw=w(l+1)/zzmass(l+1)
+             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
+          else if(w(l+1).le.0. .and. -w(l+1).le.zzmass(l)) then
+             sigw=w(l+1)/zzmass(l)
+             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
+
+!         Extended scheme (transfered mass > layer mass)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          else if(w(l+1).gt.0.) then
+             m=l+1
+             Mtot = zzmass(m)
+             MQtot = zzmass(m)*q(m)
+             do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+zzmass(m+1))))
+                m=m+1
+                Mtot = Mtot + zzmass(m)
+                MQtot = MQtot + zzmass(m)*q(m)
+             end do
+             if (m.lt.llm) then
+                sigw=(w(l+1)-Mtot)/zzmass(m+1)
+                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*(q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
+             else
+!                w(l+1) = Mtot
+!                qm(l+1) = Mqtot / Mtot
+                write(*,*) 'top layer is disappearing !',l,Mtot,w(l+1),qm(l+1)
+		print*,zzmass
+ 		print*,w
+		print*,q
+		print*,qm
+               stop
+             end if
+          else      ! if(w(l+1).lt.0) 
+             m = l-1 
+             Mtot = zzmass(m+1)
+             MQtot = zzmass(m+1)*q(m+1)
+             if (m.gt.0) then ! because some compilers will have problems
+                              ! evaluating zzmass(0)
+              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+zzmass(m))))
+                m=m-1
+                Mtot = Mtot + zzmass(m+1)
+                MQtot = MQtot + zzmass(m+1)*q(m+1)
+                if (m.eq.0) exit
+              end do
+             endif
+             if (m.gt.0) then
+                sigw=(w(l+1)+Mtot)/zzmass(m)
+                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*(q(m)-0.5*(1.+sigw)*dzq(m))  )
+             else
+                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
+             end if   
+          end if
+       enddo
+
+!     boundary conditions (not used in newcondens !!)
+!         qm(llm+1)=0.
+!         if(w(1).gt.0.) then
+!            qm(1)=q(1)
+!         else 
+!           qm(1)=0.
+!         end if
+
+      return
+      end
Index: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 728)
@@ -1,5 +1,5 @@
-subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
-
-  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD
+subroutine moistadj(pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
+
+  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
 
   implicit none
@@ -24,14 +24,14 @@
 #include "comcstfi.h"
 
-!     Pre-arguments (for universal model)
-      real pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
+      REAL pt(ngridmx,nlayermx)            ! temperature (K)
+      REAL pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
       REAL pdq(ngridmx,nlayermx,nqmx)
 
-      real dqmana(ngridmx,nlayermx,nqmx)   ! tendency of tracers (kg/kg.s-1)
-      REAL dtmana(ngridmx,nlayermx)        ! temperature increment
-
-!     Arguments
-      REAL t(ngridmx,nlayermx)       ! temperature (K)
-      REAL q(ngridmx,nlayermx)       ! humidite specifique (kg/kg)
+      REAL pdqmana(ngridmx,nlayermx,nqmx)  ! tendency of tracers (kg/kg.s-1)
+      REAL pdtmana(ngridmx,nlayermx)       ! temperature increment
+
+!     local variables
+      REAL zt(ngridmx,nlayermx)      ! temperature (K)
+      REAL zq(ngridmx,nlayermx)      ! humidite specifique (kg/kg)
       REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
       REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche (Pa)
@@ -49,5 +49,6 @@
 
 !     Local variables
-      INTEGER i, k, iq
+      INTEGER i, k, iq, nn
+      INTEGER, PARAMETER :: niter = 6
       INTEGER k1, k1p, k2, k2p
       LOGICAL itest(ngridmx)
@@ -55,8 +56,7 @@
       REAL cp_new_t(nlayermx)
       REAL cp_delta_t(nlayermx)
-      REAL new_qb(nlayermx)
       REAL v_cptj(nlayermx), v_cptjk1, v_ssig
-      REAL v_cptt(ngridmx,nlayermx), v_p, v_t
-      REAL v_qs(ngridmx,nlayermx), v_qsd(ngridmx,nlayermx)
+      REAL v_cptt(ngridmx,nlayermx), v_p, v_t, zpsat
+      REAL zqs(ngridmx,nlayermx), zdqs(ngridmx,nlayermx)
       REAL zq1(ngridmx), zq2(ngridmx)
       REAL gamcpdz(ngridmx,2:nlayermx)
@@ -94,29 +94,22 @@
 
 !     GCM -----> subroutine variables
-      DO k = 1, nlayermx
-      DO i = 1, ngridmx
-
-         q(i,k)    = pq(i,k,i_h2o)
-
-         if(q(i,k).lt.0.)then
-            q(i,k)=0.0
+      zq(1:ngridmx,1:nlayermx)    = pq(1:ngridmx,1:nlayermx,i_h2o)+ pdq(1:ngridmx,1:nlayermx,i_h2o)*ptimestep
+      zt(1:ngridmx,1:nlayermx)    = pt(1:ngridmx,1:nlayermx)
+      pdqmana(1:ngridmx,1:nlayermx,1:nqmx)=0.0
+
+      DO k = 1, nlayermx
+       DO i = 1, ngridmx
+         if(zq(i,k).lt.0.)then
+            zq(i,k)=0.0
          endif
-         DO iq = 1, nqmx
-            dqmana(i,k,iq)=0.0
-         ENDDO
-      ENDDO
-      ENDDO
-
-      DO k = 1, nlayermx
-         DO i = 1, ngridmx
-            local_q(i,k) = q(i,k)
-            local_t(i,k) = t(i,k)
-            rneb(i,k) = 0.0
-            d_ql(i,k) = 0.0
-            d_t(i,k)  = 0.0
-            d_q(i,k)  = 0.0
-         ENDDO
-         new_qb(k)=0.0
-      ENDDO
+       ENDDO
+      ENDDO
+      
+      local_q(1:ngridmx,1:nlayermx) = zq(1:ngridmx,1:nlayermx)
+      local_t(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
+      rneb(1:ngridmx,1:nlayermx) = 0.0
+      d_ql(1:ngridmx,1:nlayermx) = 0.0
+      d_t(1:ngridmx,1:nlayermx)  = 0.0
+      d_q(1:ngridmx,1:nlayermx)  = 0.0
 
 !     Calculate v_cptt
@@ -124,23 +117,11 @@
          DO i = 1, ngridmx
             v_cptt(i,k) = RCPD * local_t(i,k)
-            v_t = local_t(i,k)
+            v_t = MAX(local_t(i,k),15.)
             v_p = pplay(i,k)
 
-            call watersat(v_t,v_p,v_qs(i,k))
-            call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
+            call Psat_water(v_t,v_p,zpsat,zqs(i,k))
+	    call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
          ENDDO
       ENDDO
-
-!     TEST: RH DIAGNOSTIC
-!      DO k = 1, nlayermx
-!         DO i = 1, ngridmx
-!            v_t = local_t(i,k)
-!            IF (v_t.LT.T_h2O_ice_liq) THEN
-!               print*,'RHs=',q(i,k) / v_qs(i,k)
-!            ELSE
-!               print*,'RHl=',q(i,k) / v_qs(i,k)
-!            ENDIF
-!         ENDDO
-!      ENDDO
 
 !     Calculate Gamma * Cp * dz: (gamma is the critical gradient)
@@ -149,12 +130,8 @@
             zdp = pplev(i,k)-pplev(i,k+1)
             zdpm = pplev(i,k-1)-pplev(i,k)
-!         gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) * 
-            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) *             &
-                (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
-                +RLVTT /(zdpm+zdp) *                            &
-                (v_qs(i,k-1)*zdpm + v_qs(i,k)*zdp)              &
-                )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )     &
-                / (1.0+(v_qsd(i,k-1)*zdpm+                      &
-                v_qsd(i,k)*zdp)/(zdpm+zdp) )                    
+            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
+                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+	        * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
          ENDDO
       ENDDO
@@ -174,6 +151,6 @@
       IF (k2 .GT. nlayermx) GOTO 9999
       zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
-      zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
-         +(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
+      zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
+         +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
 
       IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
@@ -184,5 +161,5 @@
       IF (k2 .EQ. nlayermx) GOTO 821
       k2p = k2 + 1
-      zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
+      zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
       zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
 
@@ -194,9 +171,9 @@
 !------------------------------------------------------ local adjustment
   830 CONTINUE ! actual adjustment
+    Do nn=1,niter
       v_cptj(k1) = 0.0
       zdp = pplev(i,k1)-pplev(i,k1+1)
-      v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))        &
-                    + RLVTT*(local_q(i,k1)-v_qs(i,k1)) ) * zdp
-      v_ssig = zdp * (1.0+v_qsd(i,k1))
+      v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT*(local_q(i,k1)-zqs(i,k1)) ) * zdp
+      v_ssig = zdp * (1.0+zdqs(i,k1))
 
       k1p = k1 + 1
@@ -204,8 +181,6 @@
          zdp = pplev(i,k)-pplev(i,k+1)
          v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
-         v_cptjk1 = v_cptjk1 + zdp                                    &
-                  * ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k))       &
-                    + RLVTT*(local_q(i,k)-v_qs(i,k)) )        
-         v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
+         v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))*(v_cptt(i,k)+v_cptj(k)) + RLVTT*(local_q(i,k)-zqs(i,k)) )        
+         v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
       ENDDO
 
@@ -215,13 +190,19 @@
          cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
          cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
-         new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
-         local_q(i,k) = new_qb(k)
+         v_cptt(i,k)=cp_new_t(k)
+         local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
          local_t(i,k) = cp_new_t(k) / RCPD
 
-!          print*,'v_qs in loop=',v_qs
-!          print*,'v_qsd in loop=',v_qsd
-!          print*,'new_qb in loop=',new_qb
-!          print*,'cp_delta_t in loop=',cp_delta_t
-      ENDDO
+         v_t = MAX(local_t(i,k),15.)
+         v_p = pplay(i,k)
+         
+         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
+         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+
+
+
+!          print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
+      ENDDO
+    Enddo ! nn=1,niter
 
 
@@ -230,32 +211,20 @@
 !              -- the layer about to be adjusted
 
-      DO k = k1, k2
-         v_cptt(i,k) = RCPD * local_t(i,k)
-         v_t = local_t(i,k)
-         v_p = pplay(i,k)
-
-!           IF (v_t.LT.t_coup) THEN
-!              v_qs(i,k) = qsats(v_t) / v_p
-!              v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
-!           ELSE
-!              v_qs(i,k) = qsatl(v_t) / v_p
-!              v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
-!           ENDIF
-
-         call watersat(v_t,v_p,v_qs(i,k))
-         call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
-
-      ENDDO
+!      DO k = k1, k2
+!         v_cptt(i,k) = RCPD * local_t(i,k)
+!         v_t = local_t(i,k)
+!         v_p = pplay(i,k)
+!       
+!         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
+!         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+!      ENDDO
+
       DO k = 2, nlayermx
          zdpm = pplev(i,k-1) - pplev(i,k)
          zdp = pplev(i,k) - pplev(i,k+1)
-!         gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
-         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) *                       &
-                           (v_cptt(i,k-1)*zdpm+v_cptt(i,k)*zdp)        &
-                          +RLVTT /(zdpm+zdp) *                         &
-                           (v_qs(i,k-1)*zdpm+v_qs(i,k)*zdp)             &
-                         )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )    &
-                     / (1.0+(v_qsd(i,k-1)*zdpm+v_qsd(i,k)*zdp)         &
-                           /(zdpm+zdp) )
+         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)             &
+                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
       ENDDO
 
@@ -264,6 +233,6 @@
       IF (k1 .EQ. 1) GOTO 841 ! yes we have!
       zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
-      zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
-        + (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
+      zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
+        + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
       IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
 
@@ -271,5 +240,5 @@
       k1 = k1 - 1
       IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
-      zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))               &
+      zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1))               &
                   *(pplev(i,k1-1)-pplev(i,k1))
       zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
@@ -299,5 +268,5 @@
       DO i = 1, ngridmx
          IF (itest(i)) THEN
-         delta_q(i,k) = local_q(i,k) - q(i,k)
+         delta_q(i,k) = local_q(i,k) - zq(i,k)
          IF (delta_q(i,k).LT.0.) rneb(i,k)  = 1.0
          ENDIF
@@ -327,6 +296,5 @@
       DO i = 1, ngridmx
          IF (itest(i)) THEN
-         IF (zq2(i).NE.0.0) &
-           d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
+           IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
          ENDIF
       ENDDO
@@ -343,6 +311,6 @@
       DO k = 1, nlayermx
       DO i = 1, ngridmx
-         d_t(i,k) = local_t(i,k) - t(i,k)
-         d_q(i,k) = local_q(i,k) - q(i,k)
+         d_t(i,k) = local_t(i,k) - zt(i,k)
+         d_q(i,k) = local_q(i,k) - zq(i,k)
       ENDDO
       ENDDO
@@ -352,61 +320,11 @@
          DO i = 1, ngridmx
             
-            dtmana(i,k)       = d_t(i,k)/ptimestep
-            dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
-            dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
+            pdtmana(i,k)       = d_t(i,k)/ptimestep
+            pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
+            pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
          
          ENDDO
       ENDDO
 
-!      print*,'IN MANABE:'
-!      print*,'pplev=',pplev
-!      print*,'t=',t
-!      print*,'d_t=',d_t
-!      print*,'d_q=',d_q
-!      print*,'local_q=',local_q
-!      print*,'q=',q
-!      print*,'v_qs(i,k)=',v_qs
-!      print*,'v_qsd(i,k)=',v_qsd
-!      print*,'cp_delta_t(k)=',cp_delta_t
-
-!      print*,'d_ql=',d_ql
-!      print*,'delta_q=',delta_q
-!      print*,'zq1=',zq1
-!      print*,'zq2=',zq2
-!!      print*,'i_h2o=',i_h2o
-!      print*,'i_ice=',i_ice
-!
-!      print*,'IN MANABE:'
-!      print*,'d_q=',d_q
-!      print*,'d_ql=',d_ql
-!      print*,'dtmana=',d_t
-!     stop
-!      print*,'gamcpdz at end=',gamcpdz
-      !  stop    
-
-!     Some conservation diagnostics...
-!      dEtot=0.0
-!      dL1tot=0.0
-!      dL2tot=0.0
-!      dqtot=0.0
-!      masse=0.0
-!      DO k = 1, nlayermx
-!         DO i = 1, ngridmx
-!
-!            masse = (pplev(i,k) - pplev(i,k+1))/g
-!
-!            dEtot  = dEtot  + cpp*d_t(i,k)*masse
-!            dL1tot = dL1tot + RLVTT*d_ql(i,k)*masse
-!            dL2tot = dL2tot + RLVTT*d_q(i,k)*masse ! is this line necessary?
-!
-!            dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse
-!
-!         ENDDO
-!      ENDDO
-
-!        print*,'In manabe energy change=',dEtot
-!        print*,'In manabe condense energy change 1 =',dL1tot
-!        print*,'In manabe condense energy change 2 =',dL2tot
-!        print*,'In manabe water change=',dqtot
 
       RETURN
Index: trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F	(revision 728)
@@ -11,5 +11,5 @@
 #include "dimensions.h"
 #include "dimphys.h"
-!#include "comgeomfi.h"
+#include "comgeomfi.h"
 #include "surfdat.h"
 #include "planete.h"
@@ -111,53 +111,59 @@
      .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
 c
-c Read latitudes (coordinates): No need, these are provided by the dynamics
-c
-!      ierr=nf90_inq_varid(nid,"latitude",nvarid)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-!      ierr=nf90_get_var(nid,nvarid,lati)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-c
-c read longitudes (coordinates): No need, these are provided by the dynamics
-c
-!      ierr=nf90_inq_varid(nid,"longitude",nvarid)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-!      ierr=nf90_get_var(nid,nvarid,long)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-c
-c Read areas of meshes: No need, these are provided by the dynamics
-c
-!      ierr=nf90_inq_varid(nid,"area",nvarid)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Le champ <area> est absent'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-!      ierr=nf90_get_var(nid,nvarid,area)
-!      IF (ierr.NE.nf90_noerr) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         CALL abort
-!      ENDIF
-!      xmin = 1.0E+20
-!      xmax = -1.0E+20
-!      xmin = MINVAL(area)
-!      xmax = MAXVAL(area)
-!      PRINT*,'Aires des mailles <area>:', xmin, xmax
+c Lecture des latitudes (coordonnees):
+c
+      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <latitude> est absent'
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lati)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, lati)
+#endif
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
+         CALL abort
+      ENDIF
+c
+c Lecture des longitudes (coordonnees):
+c
+      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <longitude> est absent'
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, long)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, long)
+#endif
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
+         CALL abort
+      ENDIF
+c
+c Lecture des aires des mailles:
+c
+      ierr = NF_INQ_VARID (nid, "area", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <area> est absent'
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, area)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, area)
+#endif
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Lecture echouee pour <area>'
+         CALL abort
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      xmin = MINVAL(area)
+      xmax = MAXVAL(area)
+      PRINT*,'Aires des mailles <area>:', xmin, xmax
 c
 c Lecture du geopotentiel au sol:
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 728)
@@ -8,7 +8,8 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
-      use watercommon_h, only : RLVTT
+      use watercommon_h, only : RLVTT, Psat_water
       use gases_h
       use radcommon_h, only: sigma
+      use radii_mod, only: h2o_reffrad, co2_reffrad
       use aerosol_mod
       implicit none
@@ -191,5 +192,5 @@
 
       character*80 fichier 
-      integer l,ig,ierr,iq,i, tapphys,nw
+      integer l,ig,ierr,iq,iaer
 
       real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
@@ -213,6 +214,4 @@
       real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
       real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
-
-      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
 
 !     Tendencies due to various processes:
@@ -231,13 +230,20 @@
       real zdhadj(ngridmx,nlayermx)                           ! (K/s)
       real zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      real zdtmr(ngridmx,nlayermx)
       real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
       real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
       real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+      real zdumr(ngridmx,nlayermx),zdvmr(ngridmx,nlayermx)
+      real zdtsurfmr(ngridmx)
+      
+      real zdmassmr(ngridmx,nlayermx),zdpsrfmr(ngridmx)
 
       real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      real zdqevap(ngridmx,nlayermx)
       real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
       real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
       real zdqadj(ngridmx,nlayermx,nqmx)
       real zdqc(ngridmx,nlayermx,nqmx)
+      real zdqmr(ngridmx,nlayermx,nqmx),zdqsurfmr(ngridmx,nqmx)
       real zdqlscale(ngridmx,nlayermx,nqmx)
       real zdqslscale(ngridmx,nqmx)
@@ -329,6 +335,6 @@
       real qevap(ngridmx,nlayermx,nqmx)
       real tevap(ngridmx,nlayermx)
-      real dqevap(ngridmx,nlayermx)
-      real dtevap(ngridmx,nlayermx)
+      real dqevap1(ngridmx,nlayermx)
+      real dtevap1(ngridmx,nlayermx)
 
 !     included by BC for hydrology
@@ -343,5 +349,5 @@
 
 !     included by RW for RH diagnostic
-      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx)
+      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx),psat_tmp
 
 !     included by RW for hydrology
@@ -366,6 +372,6 @@
 
 !     included by BC for cloud fraction computations
-      real cloudfrac(ngridmx,nlayermx)
-      real totcloudfrac(ngridmx)
+      real,save :: cloudfrac(ngridmx,nlayermx)
+      real,save :: totcloudfrac(ngridmx)
 
 !     included by RW for vdifc water conservation test
@@ -389,18 +395,20 @@
 
 !     included by RW for variable H2O particle sizes
+      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
+      real :: reffrad1(ngridmx,nlayermx,naerkind) ! for clearsky call to callcorrk
+      real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
+      real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
+      real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
       real reffH2O(ngridmx,nlayermx)
       real reffcol(ngridmx,naerkind)
 
 !     included by RW for sourceevol
-      real ice_initial(ngridmx)!, isoil
+      real, save :: ice_initial(ngridmx)!, isoil
       real delta_ice,ice_tot
-      real ice_min(ngridmx)
-      save ice_min
-      save ice_initial
-
+      real, save :: ice_min(ngridmx)
+      
       integer num_run
-      logical ice_update
-      save ice_update
-
+      logical,save :: ice_update
+      
 !=======================================================================
       
@@ -457,7 +465,5 @@
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 
-         do ig=1,ngrid
-            albedo(ig)=albedo0(ig)
-         enddo
+         albedo(:)=albedo0(:)
 
          if(tlocked)then
@@ -473,11 +479,9 @@
          else
             print*,'WARNING! Thermal conduction in the soil turned off'
-            do ig=1,ngrid
-               capcal(ig)=1.e6
-               fluxgrd(ig)=0.
-               if(noradsurf)then
-                  fluxgrd(ig)=10.0
-               endif
-            enddo
+            capcal(:)=1.e6
+            fluxgrd(:)=0.
+            if(noradsurf)then
+                  fluxgrd(:)=10.0
+            endif
             print*,'Flux from ground = ',fluxgrd,' W m^-2'
          endif
@@ -502,7 +506,6 @@
 !        define surface as continent or ocean
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         rnat(:)=1
          do ig=1,ngridmx
-            rnat(ig)=1
-
 !            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
             if(inertiedat(ig,1).gt.1.E4)then
@@ -517,9 +520,5 @@
 !        initialise surface history variable 
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         do ig=1,ngridmx
-            do iq=1,nqmx
-               qsurf_hist(ig,iq)=qsurf(ig,iq)
-            enddo
-         enddo
+         qsurf_hist(:,:)=qsurf(:,:)
 
 !        initialise variable for dynamical heating diagnostic
@@ -556,11 +555,7 @@
             endif
 
-            do ig=1,ngrid
-!already done above (JL12)
-!               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
-               if(ice_update)then
-                  ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
-               endif
-            enddo
+            if(ice_update)then
+               ice_initial(:)=qsurf(:,igcm_h2o_ice)
+            endif
 
          endif
@@ -585,15 +580,14 @@
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      pdu(:,:) = 0.0
-      pdv(:,:) = 0.0
-!      if ( (.not.nearco2cond).and.(.not.firstcall) ) then
+      pdu(1:ngridmx,1:nlayermx) = 0.0
+      pdv(1:ngridmx,1:nlayermx) = 0.0
       if ( .not.nearco2cond ) then
-         pdt(:,:) = 0.0
-      endif ! this was the source of an evil bug...
-      pdq(:,:,:)  = 0.0
-      pdpsrf(:)   = 0.0
-      zflubid(:)  = 0.0
-      zdtsurf(:)  = 0.0
-      dqsurf(:,:) = 0.0
+         pdt(1:ngridmx,1:nlayermx) = 0.0
+      endif
+      pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+      pdpsrf(1:ngridmx)       = 0.0
+      zflubid(1:ngridmx)      = 0.0
+      zdtsurf(1:ngridmx)      = 0.0
+      dqsurf(1:ngridmx,1:nqmx)= 0.0
 
       zday=pday+ptime ! compute time, in sols (and fraction thereof)
@@ -610,13 +604,9 @@
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      do l=1,nlayer
-         do ig=1,ngrid
-            zzlay(ig,l)=pphi(ig,l)/g
-         enddo
-      enddo
-      do ig=1,ngrid
-         zzlev(ig,1)=0.
-         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
-      enddo
+      zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
+
+      zzlev(1:ngridmx,1)=0.
+      zzlev(1:ngridmx,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+
       do l=2,nlayer
          do ig=1,ngrid
@@ -631,5 +621,4 @@
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-
       do l=1,nlayer         
          do ig=1,ngrid 
@@ -657,5 +646,5 @@
               ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
 
-               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
                ztim1,ztim2,ztim3,mu0,fract)
 
@@ -707,5 +696,5 @@
                       zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
                       fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
-                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
+                      reffrad1,tau_col1,cloudfrac,totcloudfrac,             &
                       clearsky,firstcall,lastcall)
                  clearsky = .false.  ! just in case.     
@@ -722,15 +711,9 @@
                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
                     
-                    do l=1,nlayer
-                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
-                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
-                    enddo
-
-                    do nw=1,L_NSPECTV
-                       OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)		       
-                    enddo
-                    do nw=1,L_NSPECTI
-                       OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)		       
-                    enddo
+                    zdtlw(ig,1:nlayermx) = ntf*zdtlw1(ig,1:nlayermx) + tf*zdtlw(ig,1:nlayermx)
+                    zdtsw(ig,1:nlayermx) = ntf*zdtsw1(ig,1:nlayermx) + tf*zdtsw(ig,1:nlayermx)
+
+                    OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
+                    OLR_nu(ig,1:L_NSPECTV) = ntf*OLR_nu1(ig,1:L_NSPECTV) + tf*OLR_nu(ig,1:L_NSPECTV)		       
 
                  enddo
@@ -740,25 +723,14 @@
 !             Radiative flux from the sky absorbed by the surface (W.m-2)
               GSR=0.0
-              do ig=1,ngrid
-                 fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)             &
-                      +fluxsurf_sw(ig)*(1.-albedo(ig))
-
-                 if(noradsurf)then ! no lower surface; SW flux just disappears
-                    GSR = GSR + fluxsurf_sw(ig)*area(ig)
-                    fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
-                 endif
-
-              enddo
-              if(noradsurf)then
-                 print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
+              fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)*(1.-albedo(1:ngridmx))
+
+              if(noradsurf)then ! no lower surface; SW flux just disappears
+                 GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
+                 fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
+                 print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
               endif
 
 !             Net atmospheric radiative heating rate (K.s-1)
-              do l=1,nlayer
-                 do ig=1,ngrid
-                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
-                 enddo
-              enddo
-
+              dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
 
            elseif(newtonian)then
@@ -768,7 +740,5 @@
               call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
 
-              do ig=1,ngrid
-                 zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
-              enddo
+              zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
               ! e.g. surface becomes proxy for 1st atmospheric layer ?
 
@@ -777,19 +747,16 @@
 !          c) Atmosphere has no radiative effect
 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              do ig=1,ngrid
-                 fluxtop_dn(ig)  = fract(ig)*mu0(ig)*Fat1AU/dist_star**2
-                 if(ngrid.eq.1)then ! / by 4 globally in 1D case...
-                    fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
-                 endif
-                 fluxsurf_sw(ig) = fluxtop_dn(ig)
-                 fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
-                 fluxtop_lw(ig)  = emis(ig)*sigma*tsurf(ig)**4
-              enddo ! radiation skips the atmosphere entirely
-
-              do l=1,nlayer
-                 do ig=1,ngrid
-                    dtrad(ig,l)=0.0
-                 enddo
-              enddo ! hence no atmospheric radiative heating
+              fluxtop_dn(1:ngridmx)  = fract(1:ngridmx)*mu0(1:ngridmx)*Fat1AU/dist_star**2
+              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
+                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
+              endif
+              fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
+              fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
+              fluxtop_lw(1:ngridmx)  = emis(1:ngridmx)*sigma*tsurf(1:ngridmx)**4
+              ! radiation skips the atmosphere entirely
+
+
+              dtrad(1:ngridmx,1:nlayermx)=0.0
+              ! hence no atmospheric radiative heating
 
            endif                ! if corrk
@@ -801,15 +768,8 @@
 !    ------------------------------------------
 
-        do ig=1,ngrid
-           zplanck(ig)=tsurf(ig)*tsurf(ig)
-           zplanck(ig)=emis(ig)*sigma*zplanck(ig)*zplanck(ig)
-           fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
-        enddo
-
-        do l=1,nlayer
-           do ig=1,ngrid
-              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
-           enddo
-        enddo
+        zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
+        zplanck(1:ngridmx)=emis(1:ngridmx)*sigma*zplanck(1:ngridmx)*zplanck(1:ngridmx)
+        fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
+        pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
 
 !-------------------------
@@ -840,10 +800,8 @@
       if (calldifv) then
       
-         do ig=1,ngrid
-            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
-         enddo
-
-         zdum1(:,:)=0.0
-         zdum2(:,:)=0.0
+         zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
+
+         zdum1(1:ngridmx,1:nlayermx)=0.0
+         zdum2(1:ngridmx,1:nlayermx)=0.0
 
 
@@ -857,13 +815,9 @@
               zdum1,zdum2,pdt,pdq,zflubid,           &
               zdudif,zdvdif,zdtdif,zdtsdif,          &
-	      sensibFlux,q2,zdqdif,zdqsdif,lastcall)
+	      sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
 
 	 else
 	 
-           do l=1,nlayer
-             do ig=1,ngrid 
-               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
-             enddo
-           enddo
+           zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
  
            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
@@ -875,38 +829,16 @@
 	      sensibFlux,q2,zdqdif,zdqsdif,lastcall)
 
-           do l=1,nlayer
-             do ig=1,ngrid
-                zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
-             enddo
-           enddo
+           zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
+	   zdqevap(1:ngridmx,1:nlayermx)=0.
 
          end if !turbdiff
 
-         do l=1,nlayer
-            do ig=1,ngrid
-               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
-               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
-               pdt(ig,l)=pdt(ig,l)+zdtdif(ig,l)
-            enddo
-         enddo
-
-         do ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-         enddo
-
+         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
+         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
+         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
+         zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
          if (tracer) then 
-           do iq=1, nq
-            do l=1,nlayer
-              do ig=1,ngrid
-                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
-              enddo
-            enddo
-           enddo
-           do iq=1, nq
-              do ig=1,ngrid
-                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
-              enddo
-           enddo
-
+           pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx) 
+           dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
          end if ! of if (tracer)
 
@@ -965,7 +897,5 @@
          if(.not.newtonian)then
 
-            do ig=1,ngrid
-               zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
-            enddo
+            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
 
          endif
@@ -980,12 +910,8 @@
       if(calladj) then
 
-         do l=1,nlayer
-            do ig=1,ngrid
-               zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
-            enddo
-         enddo
-         zduadj(:,:)=0.0
-         zdvadj(:,:)=0.0
-         zdhadj(:,:)=0.0
+         zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
+         zduadj(1:ngridmx,1:nlayermx)=0.0
+         zdvadj(1:ngridmx,1:nlayermx)=0.0
+         zdhadj(1:ngridmx,1:nlayermx)=0.0
 
 
@@ -997,21 +923,11 @@
               zdqadj)
 
-         do l=1,nlayer
-            do ig=1,ngrid
-               pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
-               pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
-               pdt(ig,l)    = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
-               zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
-            enddo
-         enddo
-
+         pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
+         pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
+         pdt(1:ngridmx,1:nlayermx)    = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
+         zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
+ 
          if(tracer) then 
-           do iq=1, nq
-            do l=1,nlayer
-              do ig=1,ngrid
-                 pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq) 
-              enddo
-            enddo
-           enddo
+            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx) 
          end if
 
@@ -1061,22 +977,11 @@
               zdqc,reffrad)
 
-         do l=1,nlayer
-           do ig=1,ngrid
-             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
-             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
-             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
-           enddo
-         enddo
-         do ig=1,ngrid
-            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
-         enddo
-
-         do iq=1,nq ! should use new notation here !
-            do l=1,nlayer
-               do ig=1,ngrid
-                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
-               enddo
-            enddo
-         enddo
+
+         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
+         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
+         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
+         zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
+
+         pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx) 
          ! Note: we do not add surface co2ice tendency
          ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
@@ -1108,35 +1013,44 @@
 !        Water ice condensation in the atmosphere
 !        ----------------------------------------
-            if(watercond)then
-
-	       if(RLVTT.gt.1.e-8)then
-
+            if(watercond.and.(RLVTT.gt.1.e-8))then
+
+!             ----------------
+!             Moist convection
+!             ----------------
                ! Re-evaporate cloud water/ice
-               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
-               DO l = 1, nlayer  
-                  DO ig = 1, ngrid 
-                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
-                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
-                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
-                  enddo
-               enddo
-
-               call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
-               do l=1,nlayer
-                  do ig=1,ngrid
-                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
-                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
-                     pdt(ig,l)              = pdt(ig,l)+dtmoist(ig,l)
-                  enddo
-               enddo
+!	       dqevap1(1:ngridmx,1:nlayermx)=0.
+!	       dtevap1(1:ngridmx,1:nlayermx)=0.
+!               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
+               
+	       dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
+	       dtmoist(1:ngridmx,1:nlayermx)=0.
+
+               call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
+
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)   &
+	                  +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)     &
+                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
+
+!               do l=1,nlayermx
+!	         do ig=1,ngridmx
+!	          if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
+!	          if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
+!	          if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l		  
+!                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
+!	         enddo
+!	       enddo
+!	       print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
 
                !-------------------------
                ! test energy conservation
                if(enertest)then
-                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
+                  dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
                   do ig = 1, ngrid
-                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
+                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
                   enddo
                   print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
+                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
 		  
 		! test energy conservation
@@ -1146,48 +1060,33 @@
                endif
                !-------------------------
-               
-
-	       endif
 	       
 
-               ! Re-evaporate cloud water/ice
-               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
-               do l = 1, nlayer  
-                  do ig = 1, ngrid
-                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
-                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
-                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
-                  enddo
-               enddo ! note: we use qevap but not tevap in largescale/moistadj
-                     ! otherwise is a big mess
-
-               call largescale(ptimestep,pplev,pplay,pt,qevap,        & ! a bug was here!
-                    pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
-               do l=1,nlayer
-                  do ig=1,ngrid
-                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
-                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
-                     pdt(ig,l)              = pdt(ig,l)+dtlscale(ig,l)
-                     
-                     if(aeroh2o.and.(.not.aerofixh2o)) then
-                        reffrad(ig,l,iaero_h2o)=reffH2O(ig,l)
-                     endif
-
-                  enddo
-               enddo
+!        --------------------------------
+!        Large scale condensation/evaporation
+!        --------------------------------
+
+               call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
+
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
 
                !-------------------------
                ! test energy conservation
                if(enertest)then
-                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
                   do ig = 1, ngrid
-                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
+                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
                   enddo
-                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
-		  
-		! test energy conservation
-	          dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
-                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
-                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
+                  dEtot=cpp*SUM(massarea(:,:)*(dtlscale(:,:)))/totarea
+		  if(isnan(dEtot)) then
+		     print*,'Nan in largescale, abort'
+                     STOP
+		  endif
+		  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
+
+               ! test water conservation
+	          dWtot = SUM(massarea(:,:)*dqvaplscale(:,:))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea
+                  print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
                endif
                !-------------------------
@@ -1200,8 +1099,5 @@
                enddo
 
-               ! compute total cloud fraction in column
-               call totalcloudfrac(cloudfrac,totcloudfrac)
-
-           endif                ! of if (watercondense)
+            endif                ! of if (watercondense)
            
 
@@ -1209,27 +1105,21 @@
 !        Water ice / liquid precipitation
 !        --------------------------------
-           if(waterrain)then
-
-              zdqrain(:,:,:) = 0.0
-              zdqsrain(:)    = 0.0
-              zdqssnow(:)    = 0.0
-
-              call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
+            if(waterrain)then
+
+               zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+               zdqsrain(1:ngridmx)    = 0.0
+               zdqssnow(1:ngridmx)    = 0.0
+
+               call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
                    zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
 
-              do l=1,nlayer
-                 do ig=1,ngrid
-                    pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
-                    pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
-                    pdt(ig,l)              = pdt(ig,l)+zdtrain(ig,l)
-                 enddo
-              enddo
-
-              do ig=1,ngrid
-                 dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
-                 dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
-                 rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
-              enddo
-
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
+	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
+	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
+               dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
+               dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
+               rainout(1:ngridmx)             = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic   
 
 
@@ -1242,5 +1132,5 @@
 		  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
                   dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
-		  dVtot = dItot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
+		  dVtot = dVtot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
 		  dEtot = dItot + dVtot
                  print*,'In rain dItot =',dItot,' W m-2'
@@ -1258,6 +1148,6 @@
               !-------------------------
 
-           end if                 ! of if (waterrain)
-        end if                    ! of if (water)
+            end if                 ! of if (waterrain)
+         end if                    ! of if (water)
 
 
@@ -1268,6 +1158,6 @@
 !        ------------- 
         if (sedimentation) then 
-           zdqsed(:,:,:) = 0.0
-           zdqssed(:,:)  = 0.0
+           zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+           zdqssed(1:ngridmx,1:nqmx)  = 0.0
 
 
@@ -1297,15 +1187,9 @@
            !-------------------------
 
-           do iq=1,nq
-              ! for now, we only allow H2O ice to sediment
+           ! for now, we only allow H2O ice to sediment
               ! and as in rain.F90, whether it falls as rain or snow depends
               ! only on the surface temperature
-              do ig=1,ngrid
-                 do l=1,nlayer
-                    pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
-                 enddo
-                 dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
-              enddo
-           enddo
+           pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
+           dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
 
            !-------------------------
@@ -1326,4 +1210,39 @@
 !     ---------
 
+!       -----------------------------------
+!       Updating atm mass and tracer budget
+!       -----------------------------------
+
+         if(mass_redistrib) then
+
+            zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
+	       (   zdqevap(1:ngridmx,1:nlayermx)                          &
+	         + zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
+		 + dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
+		 + dqvaplscale(1:ngridmx,1:nlayermx) )
+		 
+            
+	    call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
+	    call writediagfi(ngridmx,"mass","mass"," ",3,mass)
+
+	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
+                       rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
+		       pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
+		       zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
+	 
+
+            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
+            dqsurf(1:ngridmx,1:nqmx)         = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
+            pdt(1:ngridmx,1:nlayermx)        = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
+            pdu(1:ngridmx,1:nlayermx)        = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
+            pdv(1:ngridmx,1:nlayermx)        = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
+	    pdpsrf(1:ngridmx)                = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
+            zdtsurf(1:ngridmx)               = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
+	    
+!	    print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap) 
+	 endif
+
+
+
 !       ---------------------------------
 !       Updating tracer budget on surface
@@ -1336,10 +1255,6 @@
             ! note: for now, also changes albedo in the subroutine
 
-            do ig=1,ngrid
-               zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
-               do iq=1,nq
-                  qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
-               enddo
-            enddo
+            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
+            qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
             ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
 
@@ -1365,9 +1280,5 @@
          ELSE     ! of if (hydrology)
 
-            do iq=1,nq
-               do ig=1,ngrid
-                  qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
-               enddo
-            enddo           
+            qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
 
          END IF   ! of if (hydrology)
@@ -1380,7 +1291,5 @@
 
          if(ice_update)then
-            do ig = 1, ngrid
-               ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
-            enddo
+            ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
          endif
 
@@ -1393,7 +1302,5 @@
 
 !     Increment surface temperature
-      do ig=1,ngrid
-         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
-      enddo
+      tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx) 
 
 !     Compute soil temperatures and subsurface heat flux
@@ -1406,9 +1313,5 @@
 ! test energy conservation
       if(enertest)then
-         dEtots=0.0
-         do ig = 1, ngrid
-            dEtots = dEtots + capcal(ig)*zdtsurf(ig)*area(ig)
-         enddo
-         dEtots=dEtots/totarea
+         dEtots = SUM(area(:)*capcal(:)*zdtsurf(:))/totarea
          print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
@@ -1425,49 +1328,31 @@
  
       ! temperature, zonal and meridional wind
-      do l=1,nlayer
-        do ig=1,ngrid
-          zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
-          zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
-          zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
-
-          ! diagnostic
-          zdtdyn(ig,l)     = ztprevious(ig,l)-pt(ig,l)
-          ztprevious(ig,l) = zt(ig,l)
-        enddo
+      zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
+      zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
+      zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
+
+      ! diagnostic
+      zdtdyn(1:ngridmx,1:nlayermx)     = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
+      ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
+
+      if(firstcall)then
+         zdtdyn(1:ngridmx,1:nlayermx)=0.0
+      endif
+
+      ! dynamical heating diagnostic
+      do ig=1,ngrid
+         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
       enddo
 
-      if(firstcall)then
-         zdtdyn(:,:)=0.0
-      endif
-
-      ! dynamical heating diagnostic
-      fluxdyn(:)=0.
-      do ig=1,ngrid
-         do l=1,nlayer
-            fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)   &
-                 *(pplev(ig,l)-pplev(ig,l+1))*cpp/g
-         enddo
-      enddo
-
       ! tracers
-      do iq=1, nq
-        do l=1,nlayer
-          do ig=1,ngrid
-            zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-          enddo
-        enddo
-      enddo
+      zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
 
       ! surface pressure
-      do ig=1,ngrid
-          ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
-      enddo
+      ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
 
       ! pressure
       do l=1,nlayer
-        do ig=1,ngrid
-             zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
-             zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
-        enddo
+          zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
+          zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
       enddo
 
@@ -1558,30 +1443,25 @@
 !     ---------------------------------------------------------
       if(tracer)then
-         qcol(:,:)=0.0
+         qcol(1:ngridmx,1:nqmx)=0.0
          do iq=1,nq
-            do ig=1,ngrid
-               do l=1,nlayer
-                  qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *    &
-                            (pplev(ig,l) - pplev(ig,l+1)) / g
-               enddo
+           do ig=1,ngridmx
+              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
+           enddo
+         enddo
+
+         ! Generalised for arbitrary aerosols now. (LK)
+         reffcol(1:ngridmx,1:naerkind)=0.0
+         if(co2cond.and.(iaero_co2.ne.0))then
+            call co2_reffrad(zq,reffrad)
+            do ig=1,ngridmx
+               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
             enddo
-         enddo
-
-         ! Generalised for arbitrary aerosols now. (LK)
-         reffcol(:,:)=0.0
-         do ig=1,ngrid
-            do l=1,nlayer
-               if(co2cond.and.(iaero_co2.ne.0)) then
-                  reffcol(ig,iaero_co2) = reffcol(ig,iaero_co2) + zq(ig,l,igcm_co2_ice) * &
-                                  reffrad(ig,l,iaero_co2) *    &
-                                  (pplev(ig,l) - pplev(ig,l+1)) / g
-               endif
-               if(water.and.(iaero_h2o.ne.0)) then
-                  reffcol(ig,iaero_h2o) = reffcol(ig,iaero_h2o) + zq(ig,l,igcm_h2o_ice) * &
-                                  reffrad(ig,l,iaero_h2o) *    &
-                                  (pplev(ig,l) - pplev(ig,l+1)) / g
-               endif
+         endif
+         if(water.and.(iaero_h2o.ne.0))then
+            call h2o_reffrad(zq,zt,reffrad,nueffrad)
+            do ig=1,ngridmx
+               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
             enddo
-         enddo
+         endif
 
       endif
@@ -1622,5 +1502,7 @@
          do l = 1, nlayer
             do ig = 1, ngrid
-               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
+!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
+               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
+
                RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
             enddo
@@ -1629,9 +1511,5 @@
          ! compute maximum possible H2O column amount (100% saturation)
          do ig=1,ngrid
-            H2Omaxcol(ig)=0.0
-            do l=1,nlayer
-               H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) *    &
-                    (pplev(ig,l) - pplev(ig,l+1))/g
-            enddo
+               H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
          enddo
 
@@ -1662,12 +1540,9 @@
 
                   ! add multiple years of evolution
-                  qsurf_hist(ig,igcm_h2o_ice) = &
-                     !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0 
-                     qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
+                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
 
                   ! if ice has gone -ve, set to zero
                   if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
                      qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                     !qsurf_hist(ig,igcm_h2o_vap) = 0.0
                   endif
 
@@ -1675,5 +1550,4 @@
                   if(ice_min(ig).lt.0.01)then
                      qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                     !qsurf_hist(ig,igcm_h2o_vap) = 0.0
                   endif
 
@@ -1681,11 +1555,6 @@
 
                ! enforce ice conservation
-               ice_tot=0.0
-               do ig = 1, ngrid
-                  ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
-               enddo
-               do ig = 1, ngrid
-                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
-               enddo
+               ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
+               qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
 
             endif
@@ -1803,4 +1672,5 @@
 	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
 	     call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
+	     call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
 	     call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
           endif 
@@ -1818,16 +1688,12 @@
 
 !     Output aerosols
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
-           call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3, &
-           reffrad(1,1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
-           call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3, &
-           reffrad(1,1,iaero_h2o))
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
-           call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol', &
-           'um kg m^-2',2,reffcol(1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
-           call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol', &
-           'um kg m^-2',2,reffcol(1,iaero_h2o))
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	  call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	  call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	  call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	  call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
 
 !     Output tracers
@@ -1843,7 +1709,7 @@
 
           if(watercond.or.CLFvarying)then
-!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
              call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
           endif
@@ -1877,6 +1743,6 @@
 !      output spectrum 
       if(specOLR.and.corrk)then      
-                 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
-                 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
+         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
+         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
       endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 728)
+++ trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 728)
@@ -0,0 +1,325 @@
+!==================================================================
+module radii_mod
+!==================================================================
+!  module to centralize the radii calculations for aerosols
+! OK for water but should be extended to other aerosols (CO2,...)
+!==================================================================
+      
+!     water cloud optical properties
+      
+      real, save ::  rad_h2o
+      real, save ::  rad_h2o_ice
+      real, save ::  Nmix_h2o
+      real, save ::  Nmix_h2o_ice
+      real, parameter ::  coef_chaud=0.13
+      real, parameter ::  coef_froid=0.09
+
+
+      contains
+
+
+!==================================================================
+   subroutine su_aer_radii(reffrad,nueffrad)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of liquid and icy water particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+ ! to use  'getin'
+      use ioipsl_getincom 
+      use radinc_h, only: naerkind
+      use aerosol_mod
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+
+      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind)     !variance     
+
+      logical, save :: firstcall=.true.
+      integer :: iaer   
+      
+      print*,'enter su_aer_radii'
+          do iaer=1,naerkind
+!     these values will change once the microphysics gets to work
+!     UNLESS tracer=.false., in which case we should be working with
+!     a fixed aerosol layer, and be able to define reffrad in a 
+!     .def file. To be improved!
+
+            if(iaer.eq.iaero_co2)then ! CO2 ice
+               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-4
+               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
+            endif
+
+            if(iaer.eq.iaero_h2o)then ! H2O ice
+               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
+               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
+            endif
+
+            if(iaer.eq.iaero_dust)then ! dust
+               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
+               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
+            endif
+ 
+            if(iaer.eq.iaero_h2so4)then ! H2O ice
+               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-6
+               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
+            endif
+
+            if(iaer.gt.4)then
+               print*,'Error in callcorrk, naerkind is too high (>4).'
+               print*,'The code still needs generalisation to arbitrary'
+               print*,'aerosol kinds and number.'
+               call abort
+            endif
+
+         enddo
+
+
+         if (radfixed) then
+
+	    write(*,*)"radius of H2O water particles:"
+            rad_h2o=13. ! default value
+            call getin("rad_h2o",rad_h2o)
+            write(*,*)" rad_h2o = ",rad_h2o
+
+	    write(*,*)"radius of H2O ice particles:"
+            rad_h2o_ice=35. ! default value
+            call getin("rad_h2o_ice",rad_h2o_ice)
+            write(*,*)" rad_h2o_ice = ",rad_h2o_ice
+
+	 else
+
+            write(*,*)"Number mixing ratio of H2O water particles:"
+            Nmix_h2o=1.e6 ! default value
+            call getin("Nmix_h2o",Nmix_h2o)
+            write(*,*)" Nmix_h2o = ",Nmix_h2o
+
+            write(*,*)"Number mixing ratio of H2O ice particles:"
+            Nmix_h2o_ice=Nmix_h2o ! default value
+            call getin("Nmix_h2o_ice",Nmix_h2o_ice)
+            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
+	 endif
+
+      print*,'exit su_aer_radii'
+
+   end subroutine su_aer_radii
+!==================================================================
+
+
+!==================================================================
+   subroutine h2o_reffrad(pq,pt,reffrad,nueffrad)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of liquid and icy water particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+      use radinc_h, only: naerkind
+      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
+      use aerosol_mod, only : iaero_h2o
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "comcstfi.h"
+
+      real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
+      real, intent(in) :: pt(ngridmx,nlayermx)      !temperature (K)
+      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind) ! dispersion      
+
+      integer :: ig,l
+      real zfice ,zrad,zrad_liq,zrad_ice
+      real,external :: CBRT            
+      
+
+      if (radfixed) then
+         do l=1,nlayermx
+            do ig=1,ngridmx
+               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
+               zfice = MIN(MAX(zfice,0.0),1.0)
+               reffrad(ig,l,iaero_h2o)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
+               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
+            enddo
+         enddo
+      else
+         do l=1,nlayermx
+            do ig=1,ngridmx
+               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
+               zfice = MIN(MAX(zfice,0.0),1.0)
+	       zrad_liq  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o*pi*rhowater) )
+	       zrad_ice  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
+               zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
+	       reffrad(ig,l,iaero_h2o) = min(max(zrad,1.e-6),100.e-6)
+               enddo
+            enddo      
+      end if
+
+   end subroutine h2o_reffrad
+!==================================================================
+
+
+!==================================================================
+   subroutine h2o_cloudrad(pql,reffliq,reffice)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of liquid and icy water particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+      use watercommon_h, Only: rhowater,rhowaterice
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "comcstfi.h"
+
+      real, intent(in) :: pql(ngridmx,nlayermx) !condensed water mixing ratios (kg/kg)
+      real, intent(out) :: reffliq(ngridmx,nlayermx),reffice(ngridmx,nlayermx)     !liquid and ice water particle radii (K)
+
+      real,external :: CBRT            
+      
+
+      if (radfixed) then
+	 reffliq(1:ngridmx,1:nlayermx)= rad_h2o
+	 reffice(1:ngridmx,1:nlayermx)= rad_h2o_ice
+      else
+	 reffliq(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
+	 reffliq(1:ngridmx,1:nlayermx)  = min(max(reffliq(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
+	 reffice(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+	 reffice(1:ngridmx,1:nlayermx)  = min(max(reffice(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
+      end if
+
+   end subroutine h2o_cloudrad
+!==================================================================
+
+
+
+!==================================================================
+   subroutine co2_reffrad(pq,reffrad)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of co2 ice particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+      use radinc_h, only: naerkind
+      use aerosol_mod, only : iaero_co2
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+#include "tracer.h"
+#include "comcstfi.h"
+
+      real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
+      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
+
+      integer :: ig,l
+      real :: zrad   
+      real,external :: CBRT            
+            
+      
+
+      if (radfixed) then
+         reffrad(1:ngridmx,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
+      else
+         do l=1,nlayermx
+            do ig=1,ngridmx
+               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
+               reffrad(ig,l,iaero_co2) = min(max(zrad,1.e-6),100.e-6)
+            enddo
+         enddo      
+      end if
+
+   end subroutine co2_reffrad
+!==================================================================
+
+
+
+!==================================================================
+   subroutine dust_reffrad(reffrad)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of dust particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+      use radinc_h, only: naerkind
+      use aerosol_mod, only : iaero_dust
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
+            
+      reffrad(1:ngridmx,1:nlayermx,iaero_dust) = 2.e-6 ! dust
+
+   end subroutine dust_reffrad
+!==================================================================
+
+
+!==================================================================
+   subroutine h2so4_reffrad(reffrad)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of h2so4 particles
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+      use radinc_h, only: naerkind
+      use aerosol_mod, only : iaero_h2so4
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
+                
+      reffrad(1:ngridmx,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4
+
+   end subroutine h2so4_reffrad
+!==================================================================
+
+
+
+end module radii_mod
+!==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 728)
@@ -51,6 +51,6 @@
 !		can in princple be anything: currently it's H2O.
 !
-!     NAERKIND  The number of radiatively active species
-!               (set in scatterers.h ; built when compiling with makegcm -s #)
+!     NAERKIND  The number of radiatively active aerosol types
+!
 !     NSIZEMAX  The maximum number of aerosol particle sizes
 !
@@ -76,5 +76,5 @@
       ! equivalent temperatures are 1/NTfac of these values
       integer, parameter :: NTstar = 500
-      integer, parameter :: NTstop = 9000 ! new default for all non hot Jupiter runs
+      integer, parameter :: NTstop = 15000 ! new default for all non hot Jupiter runs
       real*8, parameter :: NTfac = 1.0D+1  
       !integer, parameter :: NTstar = 1000
Index: trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 728)
@@ -2,6 +2,8 @@
 
 
-  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2
-
+! to use  'getin'
+  use ioipsl_getincom 
+  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
+  use radii_mod, only: h2o_cloudrad
   implicit none
 
@@ -15,4 +17,5 @@
 !     -------
 !     Adapted from the LMDTERRE code by R. Wordsworth (2009)
+!     Added rain vaporization in case of T>Tsat
 !     Original author Z. X. Li (1993)
 !     
@@ -51,12 +54,17 @@
       PARAMETER (seuil_neb=0.001)
 
-      REAL,PARAMETER :: cl=2.0e-4                          ! Precipitation threshold
-      REAL,PARAMETER :: ct=1./1800.                        ! Precipitation rate
+      INTEGER,save :: precip_scheme      ! id number for precipitaion scheme
+!     for simple scheme  (precip_scheme=1)
+      REAL,SAVE :: rainthreshold                ! Precipitation threshold in simple scheme
+!     for sundquist scheme  (precip_scheme=2-3)
+      REAL,SAVE :: cloud_sat                    ! Precipitation threshold in non simple scheme
+      REAL,SAVE :: precip_timescale             ! Precipitation timescale
+!     for Boucher scheme  (precip_scheme=4)
+      REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
+      REAL,PARAMETER :: Kboucher=1.19E8
+      REAL,SAVE :: c1
 
       INTEGER ninter
       PARAMETER (ninter=5)
-
-      logical simple                    ! Use very simple Emanuel scheme?
-      parameter(simple=.true.)
 
       logical evap_prec                 ! Does the rain evaporate?
@@ -68,15 +76,7 @@
       real lconvert
 
-!     for precipitation evaporation (old scheme)
-      real eeff1
-      real eeff2
-!      parameter (eeff1=0.95)
-!      parameter (eeff2=0.98)
-      parameter (eeff1=0.5)
-      parameter (eeff2=0.8)
-
 !     Local variables
       INTEGER i, k, n
-      REAL zqs(ngridmx,nlayermx), zdelta, zcor
+      REAL zqs(ngridmx,nlayermx),Tsat(ngridmx,nlayermx), zdelta, zcor
       REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt
 
@@ -85,5 +85,7 @@
       REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
 
-      real ttemp, ptemp
+      real reffh2oliq(ngridmx,nlayermx),reffh2oice(ngridmx,nlayermx)
+  
+      real ttemp, ptemp, psat_tmp
       real tnext(ngridmx,nlayermx)
 
@@ -115,4 +117,34 @@
          PRINT*, 'in rain.F, ninter=', ninter
          PRINT*, 'in rain.F, evap_prec=', evap_prec
+
+         write(*,*) "Precipitation scheme to use?"
+         precip_scheme=1 ! default value
+         call getin("precip_scheme",precip_scheme)
+         write(*,*) " precip_scheme = ",precip_scheme
+
+         if (precip_scheme.eq.1) then
+	    write(*,*) "rainthreshold in simple scheme?"
+            rainthreshold=0. ! default value
+            call getin("rainthreshold",rainthreshold)
+            write(*,*) " rainthreshold = ",rainthreshold
+
+         else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
+	    write(*,*) "cloud water saturation level in non simple scheme?"
+            cloud_sat=2.6e-4   ! default value
+            call getin("cloud_sat",cloud_sat)
+            write(*,*) " cloud_sat = ",cloud_sat
+	    write(*,*) "precipitation timescale in non simple scheme?"
+            precip_timescale=3600.  ! default value
+            call getin("precip_timescale",precip_timescale)
+            write(*,*) " precip_timescale = ",precip_timescale
+
+         else if (precip_scheme.eq.4) then
+	    write(*,*) "multiplicative constant in Boucher 95 precip scheme"
+            Cboucher=1.   ! default value
+            call getin("Cboucher",Cboucher)
+            write(*,*) " Cboucher = ",Cboucher	
+	    c1=1.00*1.097/rhowater*Cboucher*Kboucher  
+
+	 endif
 
          firstcall = .false.
@@ -158,5 +190,7 @@
             ttemp = zt(i,k) 
             ptemp = pplay(i,k)
-            call watersat(ttemp,ptemp,zqs(i,k))
+!            call watersat(ttemp,ptemp,zqs(i,k))
+	    call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
+	    call Tsat_water(ptemp,Tsat(i,k))
          ENDDO
       ENDDO
@@ -180,26 +214,28 @@
                IF (zrfl(i) .GT.0.) THEN
 
-                  zqev     = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here
-                  zqevt    = 2.0e-5*(1.0-q(i,k)/zqs(i,k))    &
-                       *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
-                  zqevt    = MAX (zqevt, 0.0)
-                  zqev     = MIN (zqev, zqevt)
-                  zqev     = MAX (zqev, 0.0)
-                  zrfln(i) = zrfl(i) - zqev
-		  zrfln(i) = max(zrfln(i),0.0)
-
-                  !zqev     = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
-                  !zqevt    = (zrfl(i)/l2c(i,k))*eeff2
-                  !zqev     = MIN (zqev, zqevt)
-                  !zrfln(i) = zrfl(i) - zqev*l2c(i,k)
-                  !zrfln(i)  = zrfl(i) - 1.5e-5*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i))
-                  !zrfln(i)  = min(zrfln(i),zrfl(i))
-                  ! this is what is actually written in the manual
-
-                  d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
-                  !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
-                  d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
-  
-                  zrfl(i)  = zrfln(i)
+                  if(zt(i,k).gt.Tsat(i,k))then
+!		     print*,'in rain',i,k,zrfl(i),zt(i,k),Tsat(i,k)
+		     zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT,zrfl(i))
+		     zrfl(i)=MAX(zrfl(i)-zqev,0.)
+		     d_q(i,k)=zqev/l2c(i,k)
+		     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
+!		     print*,zqev,zrfl(i),d_q(i,k),d_t(i,k)
+		  else
+                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! there is a bug here
+                     zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
+                        *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
+                     zqevt = MAX (zqevt, 0.0)
+                     zqev  = MIN (zqev, zqevt)
+                     zqev  = MAX (zqev, 0.0)
+                     zrfln(i)= zrfl(i) - zqev
+		     zrfln(i)= max(zrfln(i),0.0)
+
+                     d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
+                     !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
+                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
+                     zrfl(i)  = zrfln(i)
+		  end if
+		     
+
                ENDIF
             ENDDO
@@ -211,5 +247,5 @@
 
 
-         if(simple)then
+         if(precip_scheme.eq.1)then
 
             DO i = 1, ngridmx
@@ -234,6 +270,6 @@
             ENDDO
 
-         else
-
+         elseif (precip_scheme.ge.2) then
+          
             DO i = 1, ngridmx
                IF (rneb(i,k).GT.0.0) THEN
@@ -248,11 +284,15 @@
             ENDDO
 
+           SELECT CASE(precip_scheme)
+ !precip scheme from Sundquist 78
+ 	   CASE(2)
+
             DO n = 1, ninter
                DO i = 1, ngridmx
                   IF (rneb(i,k).GT.0.0) THEN
-                     zchau(i) = (ct*ptimestep/FLOAT(ninter)) * zoliq(i)      &
-                          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) * zfrac(i)
-
-                     ! this is the ONLY place where zneb, ct and cl are used
+                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
+
+                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale)) * zoliq(i)      &
+                          * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)**2)) * zfrac(i)
                      zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
                      zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
@@ -268,4 +308,56 @@
                ENDDO
             ENDDO            
+
+ !precip scheme modified from Sundquist 78 (in q**3)
+	   CASE(3)	    
+	    
+            DO n = 1, ninter
+               DO i = 1, ngridmx
+                  IF (rneb(i,k).GT.0.0) THEN
+                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
+
+                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale*cloud_sat**2)) * (zoliq(i)/zneb(i))**3  
+                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
+                          *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
+!                          * 0.1 * (1.0-zfrac(i))
+                     ztot(i)  = zchau(i) + zfroi(i)
+
+                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
+                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
+                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+
+                  ENDIF
+               ENDDO
+            ENDDO            
+
+ !precip scheme modified from Boucher 95
+	   CASE(4)
+
+            !recalculate liquid water particle radii
+	    call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
+
+            DO n = 1, ninter
+               DO i = 1, ngridmx
+                  IF (rneb(i,k).GT.0.0) THEN
+                     ! this is the ONLY place where zneb and c1 are used
+
+                     zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) &
+                                    *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 
+                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
+                          *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
+!                          * 0.1 * (1.0-zfrac(i))
+                     ztot(i)  = zchau(i) + zfroi(i)
+
+                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
+                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
+                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+
+                  ENDIF
+               ENDDO
+            ENDDO            
+
+           END SELECT ! precip_scheme 
 
             ! Change in cloud density and surface H2O values
@@ -277,5 +369,6 @@
             ENDDO
 
-         endif ! if simple
+
+         endif ! if precip_scheme=1
 
  9999 continue
Index: trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 728)
@@ -26,4 +26,5 @@
 !     B. Charnay
 !     A. Spiga
+!     J. Leconte (2012)
 !
 !==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90	(revision 728)
@@ -1,7 +1,10 @@
       subroutine su_watercycle
 
+! to use  'getin'
+      use ioipsl_getincom 
       use watercommon_h
       implicit none
 #include "comcstfi.h"
+#include "callkeys.h"
 
 
@@ -10,9 +13,10 @@
 !     Purpose
 !     -------
-!     Set up relevant constants and parameters for the water cycle
+!     Set up relevant constants and parameters for the water cycle, and water cloud properties
 !
 !     Authors
 !     -------
 !     Robin Wordsworth (2010)
+!     Jeremy Leconte (2012)
 !
 !==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 728)
@@ -1,4 +1,5 @@
-      subroutine totalcloudfrac(rneb,totalrneb)
+      subroutine totalcloudfrac(rneb,totalrneb,pplev,pq,tau)
 
+      use watercommon_h
       implicit none
 
@@ -22,10 +23,17 @@
 #include "comgeomfi.h"
 #include "comdiurn.h"
+#include "callkeys.h"
+
 
       real rneb(ngridmx,nlayermx)             ! cloud fraction     
+      real pplev(ngridmx,nlayermx+1)
+      real pq(ngridmx,nlayermx,nqmx)
+      real tau(ngridmx,nlayermx)
+
       real totalrneb(ngridmx)                 ! total cloud fraction 
-
-      integer recovery
-      parameter(recovery=1)
+      real, dimension(nlayermx+1) :: masse
+      integer, parameter          :: recovery=7
+      integer ltau_max
+      real massetot
 
 ! hypothesis behind recovery. value: 
@@ -33,9 +41,12 @@
 ! 2 = maximal recovery
 ! 3 = minimal recovery
-
-
+! 4 = fixed recovery
+! 5 = recovery on the thicker layer
 !     Local variables
       integer ig, l
-      real clear
+      real clear,tau_min
+      real, parameter ::  tau_c=0.1 !threshold of optical depth for the calculation of total cloud fraction 
+      real rneb2(nlayermx)
+
 
       do ig=1,ngridmx
@@ -51,5 +62,5 @@
          elseif (recovery.eq.2) then
             totalrneb(ig) = rneb(ig,1)
-            do l=2,nlayermx       
+            do l=2,14 !nlayermx       
                totalrneb(ig) = max(rneb(ig,l),totalrneb(ig))
             enddo
@@ -60,6 +71,38 @@
                totalrneb(ig) = min(rneb(ig,l),totalrneb(ig))
             enddo
-         endif                  ! (recovery=)
          
+         elseif (recovery.eq.4) then
+            totalrneb(ig) = CLFfixval
+
+         elseif (recovery.eq.5) then
+            totalrneb(ig) = rneb(ig,1)            
+            do l=1,nlayermx
+               masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
+            enddo
+            ltau_max=maxloc(masse,dim=1)
+            totalrneb(ig) = rneb(ig,ltau_max)
+
+         elseif (recovery.eq.6) then
+            totalrneb(ig) = 0.            
+            do l=1,nlayermx
+               masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
+               masse(l)=max(masse(l),0.)
+            enddo
+            massetot=sum(masse,dim=1)
+            do l=1,nlayermx
+               totalrneb(ig) = totalrneb(ig)+rneb(ig,l)*masse(l)/massetot
+            enddo
+
+         elseif (recovery.eq.7) then
+
+            rneb2(:)=rneb(ig,1:nlayermx)
+	    tau_min=MIN(tau_c,MAXVAL(tau(ig,1:nlayermx))/2.)
+            do l=1,nlayermx
+               if(tau(ig,l)<tau_min) rneb2(l)=0.      	
+            enddo
+            totalrneb(ig)=maxval(rneb2(1:nlayermx))
+
+         endif                  ! (recovery=)   
+
          totalrneb(ig) = min(1.,totalrneb(ig))
          totalrneb(ig) = max(0.,totalrneb(ig))
Index: trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 728)
@@ -5,7 +5,7 @@
           pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
-          pdqdif,pdqsdif,lastcall)
-
-      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
+          pdqdif,pdqevap,pdqsdif,lastcall)
+
+      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
       use radcommon_h, only : sigma
 
@@ -112,6 +112,9 @@
       INTEGER iq
       REAL zq(ngridmx,nlayermx,nqmx)
+      REAL zqnoevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
+      REAL pdqevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
+      REAL zdmassevap(ngridmx)
       REAL rho(ngridmx)         ! near-surface air density
-      REAL qsat(ngridmx)
+      REAL qsat(ngridmx),psat(ngridmx)
       DATA firstcall/.true./
       REAL kmixmin
@@ -119,8 +122,7 @@
 !     Variables added for implicit latent heat inclusion
 !     --------------------------------------------------
-      real dqsat(ngridmx), qsat_temp1, qsat_temp2
-
-      integer ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
-      save ivap, iliq, iliq_surf,iice_surf
+      real dqsat(ngridmx), qsat_temp1, qsat_temp2,psat_temp
+
+      integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
 
       real, parameter :: karman=0.4
@@ -203,4 +205,11 @@
             ENDDO
          ENDDO
+         if (water) then
+            DO ilev=1,nlay
+               DO ig=1,ngrid
+                  zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
+               ENDDO
+            ENDDO
+         Endif
       end if
 
@@ -246,5 +255,5 @@
       end if
 
-!JL12 change zkv at the surface by zcdv to calculate the surface momentum properly
+!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
       DO ig=1,ngrid
          zkv(ig,1)=zcdv(ig)
@@ -253,5 +262,5 @@
  
 !JL12 calculate the flux coefficients (tables multiplied elements by elements)
-      zfluxv=zkv(:,1:nlay)*zb0
+      zfluxv(1:ngridmx,1:nlayermx)=zkv(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx)
       
 !-----------------------------------------------------------------------
@@ -344,11 +353,11 @@
          zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
          zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
-         zkh(ig,1)=zcdh(ig)
+         zkh(ig,1)= zcdh(ig)
       ENDDO
 
 !     JL12 calculate the flux coefficients (tables multiplied elements by elements)
 !     ---------------------------------------------------------------
-      zfluxq=zkh(:,1:nlay)*zb0 !JL12 we save zfluxq which doesn't need the Exner factor
-      zfluxt=zfluxq*zExner
+      zfluxq(1:ngridmx,1:nlayermx)=zkh(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
+      zfluxt(1:ngridmx,1:nlayermx)=zfluxq(1:ngridmx,1:nlayermx)*zExner(1:ngridmx,1:nlayermx)
 
 
@@ -502,12 +511,20 @@
 
                 ! Calculate the value of qsat at the surface (water)
-                call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
-                call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
-                call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
+                call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
+                call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
+                call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
                 dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
                 ! calculate dQsat / dT by finite differences
                 ! we cannot use the updated temperature value yet...
- 
+!		if(qsat2(ig).gt.1.) then
+!		   qsat(ig)=1.
+!		   dqsat(ig)=0.
+!		else
+!		   qsat(ig)=qsat2(ig)
+!		endif
                enddo
+
+               call writediagfi(ngrid,'qsatused','saturation mixing ratio at surface','',2,qsat)
+               call writediagfi(ngrid,'psat','saturation pressure at surface','',2,psat)
 
 ! coefficients for q
@@ -632,5 +649,42 @@
            endif                ! if (water et iq=ivap)
         end do                  ! of do iq=1,nq
-      endif                     ! traceur
+
+        if (water) then  ! special case where we recompute water mixing without any evaporation.
+	                 !    The difference with the first calculation then tells us where evaporated water has gone
+
+            DO ig=1,ngrid
+               z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
+               zcq(ig,nlay)=zmass(ig,nlay)*zqnoevap(ig,nlay)*z1(ig)
+               zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
+            ENDDO 
+            
+            DO ilay=nlay-1,2,-1
+               DO ig=1,ngrid
+                  z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
+                  zcq(ig,ilay)=(zmass(ig,ilay)*zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig)
+                  zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
+               ENDDO
+            ENDDO
+
+            DO ig=1,ngrid
+               z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
+               zcq(ig,1)=(zmass(ig,1)*zqnoevap(ig,1)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
+            enddo
+
+!     Starting upward calculations for simple tracer mixing (e.g., dust)
+            do ig=1,ngrid
+               zqnoevap(ig,1)=zcq(ig,1)
+            end do
+
+            do ilay=2,nlay
+               do ig=1,ngrid
+                  zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
+               end do
+            end do
+
+         endif               ! if water
+	
+	
+      endif                     ! tracer
 
 
@@ -659,8 +713,22 @@
             enddo
          enddo
+	 if (water) then
+            do ilev = 1, nlay
+               do ig=1,ngrid
+                  pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
+               enddo
+            enddo
+	    do ig=1,ngrid
+	       zdmassevap(ig)=SUM(pdqevap(ig,:)*zmass(ig,:))*ptimestep
+	    end do	    
+	 endif
       endif 
 
       if(water)then
-      call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
+         call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
+         if (tracer) then
+            call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
+            call writediagfi(ngrid,'h2oevap','evaporated water vapor mass','kg/m2',2,zdmassevap)
+	 endif
       endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 728)
@@ -1,3 +1,3 @@
-      module watercommon_h
+module watercommon_h
 
       implicit none
@@ -11,17 +11,156 @@
       real, parameter :: RLVTT = 2.257E+6 ! Latent heat of vaporization (J kg-1) 
       real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1)
-      !real, parameter :: RLVTT = 0.0 
-      !real, parameter :: RLSTT = 0.0
 
       real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) ! entails an energy sink but better description of albedo
-      !real, parameter :: RLFTT = 0.0 
       real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3)
+      real, parameter :: rhowaterice = 9.2E+2 ! mass of water (kg/m^3)
       real, parameter :: capcal_h2o_liq = 4181.3 ! specific heat capacity of liquid water J/kg/K
       real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2)
 
-!     This was the old Martian latent heat version:
-!     lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
+      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
+      
+      contains
 
-      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
+      
+!==================================================================
+      subroutine Psat_water(T,p,psat,qsat)
 
-      end module watercommon_h
+         implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg)
+!     for a given pressure (Pa) and temperature (K)
+!     Based on the Tetens formula from L.Li physical parametrization manual
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+
+!        input
+         real, intent(in) :: T, p
+  
+!        output
+         real psat,qsat
+
+! JL12 variables for tetens formula
+         real,parameter :: Pref_solid_liquid=611.14
+         real,parameter :: Trefvaporization=35.86
+         real,parameter :: Trefsublimation=7.66
+         real,parameter :: r3vaporization=17.269
+         real,parameter :: r3sublimation=21.875
+
+! checked vs. old watersat data 14/05/2012 by JL.
+
+         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
+            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation))
+         else                 ! liquid / vapour
+            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization))
+         endif
+         if(psat.gt.p) then
+            qsat=1.
+         else
+            qsat=epsi*psat/(p-(1.-epsi)*psat)
+         endif
+         return
+      end subroutine Psat_water
+
+
+
+
+!==================================================================
+      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat)
+
+         implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute L/cp*d (q_sat)/d T
+!     for a given temperature (K)! 
+!     Based on the Tetens formula from L.Li physical parametrization manual
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+
+!        input
+         real T, p, psat, qsat
+  
+!        output
+         real dqsat
+
+! JL12 variables for tetens formula
+         real,parameter :: Pref_solid_liquid=611.14
+         real,parameter :: Trefvaporization=35.86
+         real,parameter :: Trefsublimation=7.66
+         real,parameter :: r3vaporization=17.269
+         real,parameter :: r3sublimation=21.875
+
+         real :: dummy
+
+         if (psat.gt.p) then
+	    dqsat=0.
+	    return
+	 endif
+
+         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
+            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2
+         else                 ! liquid / vapour
+            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2
+         endif
+
+         dqsat=RLVTT/RCPD*qsat**2*(p/(epsi*psat))*dummy
+         return
+      end subroutine Lcpdqsat_water
+
+
+
+
+!==================================================================
+      subroutine Tsat_water(p,Tsat)
+
+         implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the saturation temperature
+!     for a given pressure (Pa) 
+!     Based on the Tetens formula from L.Li physical parametrization manual
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+
+!        input
+         real p
+  
+!        output
+         real Tsat
+
+! JL12 variables for tetens formula
+         real,parameter :: Pref_solid_liquid=611.14
+         real,parameter :: Trefvaporization=35.86
+         real,parameter :: Trefsublimation=7.66
+         real,parameter :: r3vaporization=17.269
+         real,parameter :: r3sublimation=21.875
+
+         if (p.lt.Pref_solid_liquid) then ! solid / vapour
+            Tsat =(T_h2O_ice_liq*r3sublimation- Trefsublimation*Log(p/Pref_solid_liquid))/(r3sublimation-Log(p/Pref_solid_liquid))
+         else                 ! liquid / vapour
+            Tsat =(T_h2O_ice_liq*r3vaporization- Trefvaporization*Log(p/Pref_solid_liquid))/(r3vaporization-Log(p/Pref_solid_liquid))
+         endif
+
+         return
+      end subroutine Tsat_water
+
+
+end module watercommon_h
Index: trunk/LMDZ.GENERIC/libf/phystd/watersat.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/watersat.F90	(revision 726)
+++ trunk/LMDZ.GENERIC/libf/phystd/watersat.F90	(revision 728)
@@ -29,13 +29,17 @@
 
   if (T.lt.T_h2O_ice_liq) then ! solid / vapour
-     qsat = 100.0 * epsi * 10**(2.07023 - 0.00320991             &
+     qsat = 100.0 * 10**(2.07023 - 0.00320991             &
           * T - 2484.896 / T + 3.56654 * alog10(T))
   else                 ! liquid / vapour
-     qsat = 100.0 * epsi * 10**(23.8319 - 2948.964 / T - 5.028  &
+     qsat = 100.0 * 10**(23.8319 - 2948.964 / T - 5.028  &
           * alog10(T) - 29810.16 * exp( -0.0699382 * T)  &
           + 25.21935 * exp(-2999.924/T))
   endif
-  qsat=qsat/p
-
+!  qsat=epsi*qsat/p
+  if(qsat.gt.p) then
+     qsat=1.
+  else
+     qsat=epsi*qsat/(p-(1.-epsi)*qsat)
+  endif
   return
 end subroutine watersat
