Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 853)
+++ trunk/LMDZ.GENERIC/README	(revision 858)
@@ -842,2 +842,12 @@
 - Corrected problems with allocated arrays in start2archive and newstart since the 19/09/2012 major commit
 - Modified those programs (and iniadvtrac) so that they can be used when compiling with -t 0
+
+== 20/12/2012 == EM
+- Fixed sedimentation issue: ensure in callsedim that the correct radii are 
+  provided to newsedim and also that the updated temperature and tracer
+  mixing ratios are used to compute sedimentation.
+- Updated the way aerosol radii are considered and used; routines in radii_mod
+  (h2o_reffrad, co2_reffrad, etc.) only handle a single aerosol. The idea here
+  is that these can be called from anywhere and that the caller doesn't need to
+  have the full (naerkind size) array of aerosol radii.
+- cleanup (addition of intent(..) to routine arguments) in various routines
Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 858)
@@ -5,5 +5,5 @@
        use aerosol_mod
        USE comgeomfi_h
-       USE tracer_h
+       USE tracer_h, only: noms,rho_co2,rho_ice
                   
        implicit none
@@ -46,25 +46,25 @@
 #include "comvert.h"
 
-      INTEGER ngrid,nlayer,nq
-
-      REAL pplay(ngrid,nlayer)
-      REAL pplev(ngrid,nlayer+1)
-      REAL pq(ngrid,nlayer,nq)
-      REAL aerosol(ngrid,nlayer,naerkind)
-      REAL reffrad(ngrid,nlayer,naerkind)
-      REAL QREFvis3d(ngrid,nlayermx,naerkind)
-      REAL QREFir3d(ngrid,nlayermx,naerkind)
-
-      REAL tau_col(ngrid)
-!      REAL tauref(ngrid), tau_col(ngrid)
-
-      real cloudfrac(ngrid,nlayermx)
+      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      INTEGER,INTENT(IN) :: nq     ! number of tracers
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
+      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
+      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayermx,naerkind) ! extinction coefficient in the visible
+      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayermx,naerkind)
+      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
+      ! BENJAMIN MODIFS
+      real,intent(in) :: cloudfrac(ngrid,nlayermx) ! cloud fraction
+      real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
+      logical,intent(in) :: clearsky
+
       real aerosol0
 
       INTEGER l,ig,iq,iaer
 
-      LOGICAL firstcall
-      DATA firstcall/.true./
-      SAVE firstcall
+      LOGICAL,SAVE :: firstcall=.true.
 
       REAL CBRT
@@ -79,8 +79,6 @@
       REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
       REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
-      ! BENJAMIN MODIFS
+
       real CLFtot
-      real totcloudfrac(ngrid)
-      logical clearsky
 
       ! identify tracers
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 858)
@@ -1,9 +1,9 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,emis,mu0,pplev,pplay,pt,                      & 
-          tsurf,fract,dist_star,aerosol,muvar,           &
+          tsurf,fract,dist_star,aerosol,muvar,                 &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
           OLR_nu,OSR_nu,                                       &
-          reffrad,nueffrad,tau_col,cloudfrac,totcloudfrac,     &
+          tau_col,cloudfrac,totcloudfrac,                      &
           clearsky,firstcall,lastcall)
 
@@ -15,6 +15,5 @@
       use gases_h
       use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
-      use aerosol_mod
-
+      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4
       USE tracer_h
 
@@ -45,17 +44,36 @@
 !     Layer #nlayermx is the layer at the top.
 
-!     INPUT
-      INTEGER icount
-      INTEGER ngrid,nlayer
-      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
-      REAL albedo(ngrid)                    ! SW albedo
-      REAL emis(ngrid)                      ! LW emissivity
-      REAL pplay(ngrid,nlayermx)            ! pres. level in GCM mid of layer
-      REAL pplev(ngrid,nlayermx+1)          ! pres. level at GCM layer boundaries
-
-      REAL pt(ngrid,nlayermx)               ! air temperature (K)
-      REAL tsurf(ngrid)                     ! surface temperature (K)
-      REAL dist_star,mu0(ngrid)             ! distance star-planet (AU)
-      REAL fract(ngrid)                     ! fraction of day
+      INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      integer,intent(in) :: nq ! number of tracers
+      REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg.m-2)
+      REAL,INTENT(IN) :: albedo(ngrid)   ! SW albedo
+      REAL,INTENT(IN) :: emis(ngrid)     ! LW emissivity
+      real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle
+      REAL,INTENT(IN) :: pplev(ngrid,nlayermx+1)  ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayermx)    ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pt(ngrid,nlayermx)  ! air temperature (K)
+      REAL,INTENT(IN) :: tsurf(ngrid)        ! surface temperature (K)
+      REAL,INTENT(IN) :: fract(ngrid)        ! fraction of day
+      REAL,INTENT(IN) :: dist_star           ! distance star-planet (AU)
+      REAL,INTENT(OUT) :: aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
+      real,intent(in) :: muvar(ngrid,nlayermx+1)
+      REAL,INTENT(OUT) :: dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
+      REAL,INTENT(OUT) :: dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
+      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
+      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
+      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
+      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
+      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)    ! incident top of atmosphere SW flux (W/m2)
+      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
+      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
+      REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity
+!     for H2O cloud fraction in aeropacity
+      real,intent(in) :: cloudfrac(ngrid,nlayermx)
+      real,intent(out) :: totcloudfrac(ngrid)
+      logical,intent(in) :: clearsky
+      logical,intent(in) :: firstcall ! signals first call to physics
+      logical,intent(in) :: lastcall ! signals last call to physics
 
 !     Globally varying aerosol optical properties on GCM grid
@@ -72,18 +90,6 @@
 !      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind) ! not sure of the point of these...
 
-      !!! this is a local instance of a variable saved in physiq.F and being an argument
-      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
-      REAL, INTENT(INOUT) :: nueffrad(ngrid,nlayermx,naerkind)
-
-!!     OUTPUT
-      REAL dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
-      REAL dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
-      REAL fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
-      REAL fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
-      REAL fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
-      REAL fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
-      REAL fluxtop_dn(ngrid)    ! incident top of atmosphere SW flux (W/m2)
-      REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
-      REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
+      REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
+      REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
 
 !-----------------------------------------------------------------------
@@ -117,4 +123,5 @@
 
       INTEGER ig,l,k,nw,iaer,irad
+      INTEGER icount
 
       real fluxtoplanet
@@ -126,18 +133,12 @@
 
       real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
-      REAL pq(ngrid,nlayer,nq)
-      REAL qsurf(ngrid,nq)       ! tracer on surface (e.g. kg.m-2)
-      integer nq
 
 !     Local aerosol optical properties for each column on RADIATIVE grid
-      real*8  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-
-      save qxvaer, qsvaer, gvaer
-      save qxiaer, qsiaer, giaer
+      real*8,save ::  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
+      real*8,save ::  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
+      real*8,save ::  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
+      real*8,save ::  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
+      real*8,save ::  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
+      real*8,save ::  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
 
       !REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
@@ -147,8 +148,7 @@
       real, dimension(:,:,:), save, allocatable :: QREFir3d
 
-      REAL tau_col(ngrid) ! diagnostic from aeropacity
 
 !     Misc.
-      logical firstcall, lastcall, nantest
+      logical nantest
       real*8  tempv(L_NSPECTV)
       real*8  tempi(L_NSPECTI)
@@ -169,8 +169,4 @@
       real*8 NFLUXGNDV_nu(L_NSPECTV)
 
-!     for H2O cloud fraction in aeropacity
-      real cloudfrac(ngrid,nlayermx)
-      real totcloudfrac(ngrid)
-      logical clearsky
 
       ! for weird cloud test
@@ -179,9 +175,7 @@
       real maxrad, minrad
             
-      real CBRT
-      external CBRT
+      real,external :: CBRT
 
 !     included by RW for runaway greenhouse 1D study
-      real muvar(ngrid,nlayermx+1)
       real vtmp(nlayermx)
       REAL*8 muvarrad(L_LEVELS)
@@ -208,4 +202,7 @@
 !--------------------------------------------------
 !     Effective radius and variance of the aerosols
+         allocate(reffrad(ngrid,nlayer,naerkind))
+         allocate(nueffrad(ngrid,nlayer,naerkind))
+
          if(naerkind.gt.4)then
             print*,'Code not general enough to deal with naerkind > 4 yet.'
@@ -262,7 +259,5 @@
          endif
 
-         firstcall=.false.   
-
-      end if
+      end if ! of if (firstcall)
 
 !=======================================================================
@@ -274,19 +269,20 @@
 
          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
-	    call co2_reffrad(ngrid,nq,pq,reffrad)
+	    call co2_reffrad(ngrid,nq,pq,reffrad(1,1,iaero_co2))
             print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
             print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,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(ngrid,nq,pq,pt,reffrad,nueffrad)
+	    call h2o_reffrad(ngrid,pq(1,1,igcm_h2o_ice),pt, &
+                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
             print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
          endif
          if(iaer.eq.iaero_dust)then
-	    call dust_reffrad(ngrid,reffrad)
+	    call dust_reffrad(ngrid,reffrad(1,1,iaero_dust))
             print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          if(iaer.eq.iaero_h2so4)then
-	    call h2so4_reffrad(ngrid,reffrad)
+	    call h2so4_reffrad(ngrid,reffrad(1,1,iaero_h2so4))
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
@@ -638,4 +634,6 @@
             print*,'Minimum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tgasmin=",tgasmin
             !print*,'WARNING, OVERRIDING FOR TEST'
             call abort
Index: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 858)
@@ -1,7 +1,10 @@
       SUBROUTINE callsedim(ngrid,nlay, ptimestep,
-     $                pplev,zlev, pt,
-     &                pq, pdqfi, pdqsed,pdqs_sed,nq,rfall)
+     &                pplev,zlev, pt, pdt,
+     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
 
-      USE tracer_h
+      use radinc_h, only : naerkind
+      use radii_mod, only: h2o_reffrad
+      use aerosol_mod, only : iaero_h2o
+      USE tracer_h, only : igcm_co2_ice,igcm_h2o_ice,radius,rho_q
 
       IMPLICIT NONE
@@ -33,18 +36,18 @@
 c   ----------
 
-      INTEGER ngrid              ! number of horizontal grid points
-      INTEGER nlay               ! number of atmospheric layers
-      REAL ptimestep             ! physics time step (s)
-      REAL pplev(ngrid,nlay+1)   ! pressure at inter-layers (Pa)
-      REAL pt(ngrid,nlay)        ! temperature at mid-layer (K)
-      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
-
-
-c    Traceurs :
-      integer nq             ! number of tracers
-      real pq(ngrid,nlay,nq)  ! tracers (kg/kg)
-      real pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
-      real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
-      real pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
+      integer,intent(in):: ngrid ! number of horizontal grid points
+      integer,intent(in):: nlay  ! number of atmospheric layers
+      real,intent(in):: ptimestep  ! physics time step (s)
+      real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
+      real,intent(in):: pt(ngrid,nlay)      ! temperature at mid-layer (K)
+      real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature
+      real,intent(in):: zlev(ngrid,nlay+1)  ! altitude at layer boundaries
+      integer,intent(in) :: nq ! number of tracers
+      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
+      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
+                                               ! sedimentation (kg/kg.s-1)
+      
+      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
+      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
       
 c   local:
@@ -53,15 +56,17 @@
       INTEGER l,ig, iq
 
-      real zqi(ngrid,nlayermx,nq) ! to locally store tracers
-      real masse (ngrid,nlayermx) ! Layer mass (kg.m-2)
-      real epaisseur (ngrid,nlayermx) ! Layer thickness (m)
-      real wq(ngrid,nlayermx+1) ! displaced tracer mass (kg.m-2)
+      ! for particles with varying radii:
+      real reffrad(ngrid,nlay,naerkind) ! particle radius (m) 
+      real nueffrad(ngrid,nlay,naerkind) ! aerosol effective radius variance
+
+      real zqi(ngrid,nlay,nq) ! to locally store tracers
+      real zt(ngrid,nlay) ! to locally store temperature (K)
+      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
+      real epaisseur (ngrid,nlay) ! Layer thickness (m)
+      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
 c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
-      real rfall(ngrid,nlayermx)
 
 
-      LOGICAL firstcall
-      SAVE firstcall
-      DATA firstcall/.true./
+      LOGICAL,SAVE :: firstcall=.true.
 
 c    ** un petit test de coherence
@@ -80,4 +85,5 @@
           masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
           epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
+          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
         end do
       end do
@@ -88,8 +94,7 @@
 !  of general tracers.]
  
-      zqi(1:ngrid,1:nlay,1:nq) = 0.
       do iq=1,nq
        if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
-!         (no sedim for gases; co2_ice sedim is done elsewhere)      
+!         (no sedim for gases, and co2_ice sedim is done in condense_cloud)      
 
 ! store locally updated tracers
@@ -106,12 +111,17 @@
 ! Water          
           if (water.and.(iq.eq.igcm_h2o_ice)) then
+            ! compute radii for h2o_ice 
+             call h2o_reffrad(ngrid,zqi(1,1,igcm_h2o_ice),zt,
+     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
+            ! call sedimentation for h2o_ice
              call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
-     &            pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq)
+     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
+     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq)
 
 ! General Case
           else 
              call newsedim(ngrid,nlay,1,ptimestep,
-     &            pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
-     &            zqi,wq)
+     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
+     &            zqi(1,1,iq),wq)
           endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 858)
@@ -4,5 +4,5 @@
           piceco2,psolaralb,pemisurf, &
           pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
-          pdqc,reffrad)
+          pdqc)
 
       use radinc_h, only : naerkind
@@ -65,26 +65,31 @@
 !     Arguments
 
-      INTEGER ngrid, nlayer, nq
-
-      REAL ptimestep 
-      REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
-      REAL pcapcal(ngrid)
-      REAL pt(ngrid,nlayer)
-      REAL ptsrf(ngrid)
-      REAL pphi(ngrid,nlayer)
-      REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
-      REAL pdtsrfc(ngrid),pdpsrf(ngrid)
-!      REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
-      REAL piceco2(ngrid),psolaralb(ngrid),pemisurf(ngrid)
-
-      REAL pu(ngrid,nlayer) ,  pv(ngrid,nlayer)
-      REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
-      REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
-      REAL pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
-      REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq)
-
-      REAL reffrad(ngrid,nlayer,naerkind)
-
-
+      INTEGER,INTENT(IN) :: ngrid
+      INTEGER,INTENT(IN) :: nlayer
+      INTEGER,INTENT(IN) :: nq
+      REAL,INTENT(IN) :: ptimestep 
+      REAL,INTENT(IN) :: pcapcal(ngrid)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
+      REAL,INTENT(IN) :: ptsrf(ngrid)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer)
+      REAL,INTENT(IN) :: pphi(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdt(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdu(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdv(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdtsrf(ngrid)
+      REAL,INTENT(IN) :: pu(ngrid,nlayer)
+      REAL,INTENT(IN) :: pv(ngrid,nlayer)
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
+      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
+      REAL,INTENT(INOUT) :: piceco2(ngrid)
+      REAL,INTENT(OUT) :: psolaralb(ngrid)
+      REAL,INTENT(OUT) :: pemisurf(ngrid)
+      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdtsrfc(ngrid)
+      REAL,INTENT(OUT) :: pdpsrf(ngrid)
+      REAL,INTENT(OUT) :: pduc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
 
 !-----------------------------------------------------------------------
@@ -93,4 +98,5 @@
       INTEGER l,ig,icap,ilay,it,iq
 
+      REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
       REAL*8 zt(ngrid,nlayermx)
       REAL zq(ngrid,nlayermx,nq)
@@ -123,22 +129,13 @@
 !     Saved local variables
 
-      REAL latcond
-      REAL ccond
-      REAL cpice
+      REAL,SAVE :: latcond=5.9e5
+      REAL,SAVE :: ccond
+      REAL,SAVE :: cpice=1000.
       REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
-      SAVE cpice
-      SAVE latcond,ccond
-
-      LOGICAL firstcall
-      SAVE firstcall
-      REAL SSUM
-      EXTERNAL SSUM
-
-      DATA latcond /5.9e5/
-      DATA cpice /1000./
-      DATA firstcall/.true./
-
-      REAL CBRT
-      EXTERNAL CBRT
+
+      LOGICAL,SAVE :: firstcall=.true.
+      REAL,EXTERNAL :: SSUM
+
+      REAL,EXTERNAL :: CBRT
 
       INTEGER,SAVE :: i_co2ice=0      ! co2 ice
@@ -154,8 +151,8 @@
 !     Initializations
 
-      call zerophys(ngrid*nlayer*nq,pdqc)
-      call zerophys(ngrid*nlayer*nq,pdtc)
-      call zerophys(ngrid*nlayermx*nq,zq)
-      call zerophys(ngrid*nlayermx,zt)
+      pdqc(1:ngrid,1:nlayer,1:nq)=0
+      pdtc(1:ngrid,1:nlayer)=0
+      zq(1:ngrid,1:nlayer,1:nq)=0
+      zt(1:ngrid,1:nlayer)=0
 
 !     Initialisation
@@ -320,5 +317,5 @@
             DO ig=1,ngrid
 
-               reff = reffrad(ig,ilay,iaero_co2)
+               reff = reffrad(ig,ilay)
 
                call stokes                      &
Index: trunk/LMDZ.GENERIC/libf/phystd/newsedim.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/newsedim.F	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/newsedim.F	(revision 858)
@@ -23,17 +23,17 @@
 !   ---------
 
-      INTEGER ngrid,nlay,naersize
-      REAL ptimestep            ! pas de temps physique (s)
-      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
-      REAL pt(ngrid,nlay)    ! temperature au centre des couches (K)
-      real masse (ngrid,nlay) ! masse d'une couche (kg)
-      real epaisseur (ngrid,nlay)  ! epaisseur d'une couche (m)
-      real rd(naersize)             ! particle radius (m)
-      real rho             ! particle density (kg.m-3)
-
-
-c    Traceurs :
-      real pqi(ngrid,nlay)  ! traceur   (e.g. ?/kg)
-      real wq(ngrid,nlay+1)  ! flux de traceur durant timestep (?/m-2)
+      integer,intent(in) :: ngrid  ! number of atmospheric columns
+      integer,intent(in) :: nlay  ! number of atmospheric layers
+      integer,intent(in) :: naersize   ! number of particle sizes (1 or number
+                                       ! of grid boxes)
+      real,intent(in) :: ptimestep ! physics time step (s)
+      real,intent(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressures (Pa)
+      real,intent(in) :: pt(ngrid,nlay)    ! mid-layer temperatures (K)
+      real,intent(in) :: masse (ngrid,nlay) ! atmospheric mass (kg)
+      real,intent(in) :: epaisseur (ngrid,nlay)  ! thickness of atm. layers (m)
+      real,intent(in) :: rd(naersize) ! particle radius (m)
+      real,intent(in) :: rho ! particle density (kg.m-3)
+      real,intent(inout) :: pqi(ngrid,nlay)  ! tracer   (e.g. ?/kg)
+      real,intent(out) :: wq(ngrid,nlay+1)  ! flux of tracer during timestep (?/m-2)
       
 c   local:
@@ -43,6 +43,5 @@
       REAL rfall
 
-      LOGICAL firstcall
-      SAVE firstcall
+      LOGICAL,SAVE :: firstcall=.true.
 
 c    Traceurs :
@@ -53,17 +52,14 @@
       real ptop, dztop, Ep, Stra
 
-      DATA firstcall/.true./
 
 c    Physical constant
 c    ~~~~~~~~~~~~~~~~~
-      REAL visc, molrad
 c     Gas molecular viscosity (N.s.m-2)
-      data visc/1.e-5/       ! CO2
+      real,parameter :: visc=1.e-5       ! CO2
 c     Effective gas molecular radius (m)
-      data molrad/2.2e-10/   ! CO2
+      real,parameter :: molrad=2.2e-10   ! CO2
 
 c     local and saved variable
-      real a,b
-      save a,b
+      real,save :: a,b
 
 c    ** un petit test de coherence
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 858)
@@ -12,5 +12,5 @@
       use gases_h
       use radcommon_h, only: sigma
-      use radii_mod, only: h2o_reffrad, co2_reffrad
+      use radii_mod, only: h2o_reffrad, co2_reffrad, h2o_cloudrad
       use aerosol_mod
       use surfdat_h
@@ -132,27 +132,33 @@
 !   inputs:
 !   -------
-      integer ngrid,nlayer,nq
-      real ptimestep
-      real pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      real pphi(ngrid,nlayer)
-      real pu(ngrid,nlayer),pv(ngrid,nlayer)
-      real pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
-      real pw(ngrid,nlayer)        ! pvervel transmitted by dyn3d
-      real zh(ngrid,nlayermx)      ! potential temperature (K)
-      logical firstcall,lastcall
-
-      real pday
-      real ptime 
-      logical tracerdyn
-
-      character*20 nametrac(nq)   ! name of the tracer from dynamics
+      integer,intent(in) :: ngrid ! number of atmospheric columns
+      integer,intent(in) :: nlayer ! number of atmospheric layers
+      integer,intent(in) :: nq ! number of tracers
+      character*20,intent(in) :: nametrac(nq) ! name of the tracer from dynamics
+      logical,intent(in) :: firstcall ! signals first call to physics
+      logical,intent(in) :: lastcall ! signals last call to physics
+      real,intent(in) :: pday ! number of elapsed sols since reference Ls=0
+      real,intent(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
+      real,intent(in) :: ptimestep ! physics timestep (s)
+      real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
+      real,intent(in) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
+      real,intent(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
+      real,intent(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
+      real,intent(in) :: pt(ngrid,nlayer) ! temperature (K)
+      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      real,intent(in) :: pw(ngrid,nlayer)    ! vertical velocity (m/s)
+
+
 
 !   outputs:
 !   --------
 !     physical tendencies
-      real pdu(ngrid,nlayer),pdv(ngrid,nlayer)
-      real pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
-      real pdpsrf(ngrid) ! surface pressure tendency
-
+      real,intent(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
+      real,intent(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
+      real,intent(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
+      real,intent(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
+      real,intent(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
+      logical,intent(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
 
 ! Local saved variables:
@@ -190,4 +196,5 @@
 !     for the "naerkind" optically active aerosols:
       real aerosol(ngrid,nlayermx,naerkind) 
+      real zh(ngrid,nlayermx)      ! potential temperature (K)
 
       character*80 fichier 
@@ -390,8 +397,8 @@
       real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
       real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
-      real :: nueffrad_dummy(ngrid,nlayermx,naerkind) !! AS. This is temporary. Check below why.
+!      real :: nueffrad_dummy(ngrid,nlayermx,naerkind) !! AS. This is temporary. Check below why.
       real :: reffh2oliq(ngrid,nlayermx) ! liquid water particles effective radius (m)
       real :: reffh2oice(ngrid,nlayermx) ! water ice particles effective radius (m)
-      real reffH2O(ngrid,nlayermx)
+   !   real reffH2O(ngrid,nlayermx)
       real reffcol(ngrid,naerkind)
 
@@ -406,42 +413,4 @@
 !=======================================================================
 
-      !!! ALLOCATE SAVED ARRAYS
-      IF (.not. ALLOCATED(tsurf)) ALLOCATE(tsurf(ngrid))
-      IF (.not. ALLOCATED(tsoil)) ALLOCATE(tsoil(ngrid,nsoilmx))    
-      IF (.not. ALLOCATED(albedo)) ALLOCATE(albedo(ngrid))          
-      IF (.not. ALLOCATED(albedo0)) ALLOCATE(albedo0(ngrid))          
-      IF (.not. ALLOCATED(rnat)) ALLOCATE(rnat(ngrid))         
-      IF (.not. ALLOCATED(emis)) ALLOCATE(emis(ngrid))   
-      IF (.not. ALLOCATED(dtrad)) ALLOCATE(dtrad(ngrid,nlayermx))
-      IF (.not. ALLOCATED(fluxrad_sky)) ALLOCATE(fluxrad_sky(ngrid))
-      IF (.not. ALLOCATED(fluxrad)) ALLOCATE(fluxrad(ngrid))    
-      IF (.not. ALLOCATED(capcal)) ALLOCATE(capcal(ngrid))    
-      IF (.not. ALLOCATED(fluxgrd)) ALLOCATE(fluxgrd(ngrid))  
-      IF (.not. ALLOCATED(qsurf)) ALLOCATE(qsurf(ngrid,nq))  
-      IF (.not. ALLOCATED(q2)) ALLOCATE(q2(ngrid,nlayermx+1)) 
-      IF (.not. ALLOCATED(ztprevious)) ALLOCATE(ztprevious(ngrid,nlayermx))
-      IF (.not. ALLOCATED(cloudfrac)) ALLOCATE(cloudfrac(ngrid,nlayermx))
-      IF (.not. ALLOCATED(totcloudfrac)) ALLOCATE(totcloudfrac(ngrid))
-      IF (.not. ALLOCATED(qsurf_hist)) ALLOCATE(qsurf_hist(ngrid,nq))
-      IF (.not. ALLOCATED(reffrad)) ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
-      IF (.not. ALLOCATED(nueffrad)) ALLOCATE(nueffrad(ngrid,nlayermx,naerkind))
-      IF (.not. ALLOCATED(ice_initial)) ALLOCATE(ice_initial(ngrid))
-      IF (.not. ALLOCATED(ice_min)) ALLOCATE(ice_min(ngrid)) 
-
-      IF (.not. ALLOCATED(fluxsurf_lw)) ALLOCATE(fluxsurf_lw(ngrid))
-      IF (.not. ALLOCATED(fluxsurf_sw)) ALLOCATE(fluxsurf_sw(ngrid))
-      IF (.not. ALLOCATED(fluxtop_lw)) ALLOCATE(fluxtop_lw(ngrid))
-      IF (.not. ALLOCATED(fluxabs_sw)) ALLOCATE(fluxabs_sw(ngrid))
-      IF (.not. ALLOCATED(fluxtop_dn)) ALLOCATE(fluxtop_dn(ngrid))
-      IF (.not. ALLOCATED(fluxdyn)) ALLOCATE(fluxdyn(ngrid))
-      IF (.not. ALLOCATED(OLR_nu)) ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
-      IF (.not. ALLOCATED(OSR_nu)) ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
-      IF (.not. ALLOCATED(sensibFlux)) ALLOCATE(sensibFlux(ngrid))
-      IF (.not. ALLOCATED(zdtlw)) ALLOCATE(zdtlw(ngrid,nlayermx))
-      IF (.not. ALLOCATED(zdtsw)) ALLOCATE(zdtsw(ngrid,nlayermx))
-      IF (.not. ALLOCATED(tau_col)) ALLOCATE(tau_col(ngrid))
-
-!=======================================================================
-
 ! 1. Initialisation
 ! -----------------
@@ -450,4 +419,39 @@
 !  ---------------------------------------
       if (firstcall) then
+
+        !!! ALLOCATE SAVED ARRAYS
+        ALLOCATE(tsurf(ngrid))
+        ALLOCATE(tsoil(ngrid,nsoilmx))    
+        ALLOCATE(albedo(ngrid))          
+        ALLOCATE(albedo0(ngrid))          
+        ALLOCATE(rnat(ngrid))         
+        ALLOCATE(emis(ngrid))   
+        ALLOCATE(dtrad(ngrid,nlayermx))
+        ALLOCATE(fluxrad_sky(ngrid))
+        ALLOCATE(fluxrad(ngrid))    
+        ALLOCATE(capcal(ngrid))    
+        ALLOCATE(fluxgrd(ngrid))  
+        ALLOCATE(qsurf(ngrid,nq))  
+        ALLOCATE(q2(ngrid,nlayermx+1)) 
+        ALLOCATE(ztprevious(ngrid,nlayermx))
+        ALLOCATE(cloudfrac(ngrid,nlayermx))
+        ALLOCATE(totcloudfrac(ngrid))
+        ALLOCATE(qsurf_hist(ngrid,nq))
+        ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
+        ALLOCATE(nueffrad(ngrid,nlayermx,naerkind))
+        ALLOCATE(ice_initial(ngrid))
+        ALLOCATE(ice_min(ngrid)) 
+        ALLOCATE(fluxsurf_lw(ngrid))
+        ALLOCATE(fluxsurf_sw(ngrid))
+        ALLOCATE(fluxtop_lw(ngrid))
+        ALLOCATE(fluxabs_sw(ngrid))
+        ALLOCATE(fluxtop_dn(ngrid))
+        ALLOCATE(fluxdyn(ngrid))
+        ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
+        ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
+        ALLOCATE(sensibFlux(ngrid))
+        ALLOCATE(zdtlw(ngrid,nlayermx))
+        ALLOCATE(zdtsw(ngrid,nlayermx))
+        ALLOCATE(tau_col(ngrid))
 
          !! this is defined in comsaison_h
@@ -710,5 +714,5 @@
                       zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
                       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
-                      reffrad,nueffrad,tau_col,cloudfrac,totcloudfrac,   &
+                      tau_col,cloudfrac,totcloudfrac,                    &
                       clearsky,firstcall,lastcall)      
 
@@ -724,5 +728,5 @@
                       zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
                       fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
-                      reffrad,nueffrad,tau_col1,cloudfrac,totcloudfrac,    &
+                      tau_col1,cloudfrac,totcloudfrac,                     &
                       clearsky,firstcall,lastcall)
                  clearsky = .false.  ! just in case.     
@@ -1003,5 +1007,5 @@
               qsurf(1,igcm_co2_ice),albedo,emis,            &
               zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
-              zdqc,reffrad)
+              zdqc)
 
 
@@ -1189,5 +1193,5 @@
 
            call callsedim(ngrid,nlayer,ptimestep,           &
-                pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
+                pplev,zzlev,pt,pdt,pq,pdq,zdqsed,zdqssed,nq)
 
            !-------------------------
@@ -1468,5 +1472,5 @@
          reffcol(1:ngrid,1:naerkind)=0.0
          if(co2cond.and.(iaero_co2.ne.0))then
-            call co2_reffrad(ngrid,nq,zq,reffrad)
+            call co2_reffrad(ngrid,nq,zq,reffrad(1,1,iaero_co2))
             do ig=1,ngrid
                reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
@@ -1474,11 +1478,6 @@
          endif
          if(water.and.(iaero_h2o.ne.0))then
-
-            !! AS: This modifies reffrad and nueffrad
-            !!     ... and those are saved in physiq
-            !!     To be conservative with previous versions,
-            !!     ... I put nueffrad_dummy so that nueffrad
-            !!     ... is not modified, but shouldn't it be modified actually?
-            call h2o_reffrad(ngrid,nq,zq,zt,reffrad,nueffrad_dummy)
+            call h2o_reffrad(ngrid,zq(1,1,igcm_h2o_ice),zt, &
+                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             do ig=1,ngrid
                reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
Index: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 858)
@@ -35,15 +35,15 @@
       use radinc_h, only: naerkind
       use aerosol_mod
-      USE tracer_h
-      Implicit none
-
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-
-      integer :: ngrid
-
-      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
-      real, intent(inout) :: nueffrad(ngrid,nlayermx,naerkind)     !variance     
+!      USE tracer_h
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+
+      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(out) :: nueffrad(ngrid,nlayermx,naerkind)     !variance     
 
       logical, save :: firstcall=.true.
@@ -119,5 +119,6 @@
 
 !==================================================================
-   subroutine h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
+!   subroutine h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
+   subroutine h2o_reffrad(ngrid,pq,pt,reffrad,nueffrad)
 !==================================================================
 !     Purpose
@@ -130,8 +131,8 @@
 !
 !==================================================================
-      use radinc_h, only: naerkind
+!      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
-      USE tracer_h
+!      use aerosol_mod, only : iaero_h2o
+!      USE tracer_h, only: igcm_h2o_ice
       Implicit none
 
@@ -141,10 +142,14 @@
 #include "comcstfi.h"
 
-      integer :: ngrid,nq
-
-      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
-      real, intent(in) :: pt(ngrid,nlayermx)      !temperature (K)
-      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
-      real, intent(inout) :: nueffrad(ngrid,nlayermx,naerkind) ! dispersion      
+      integer,intent(in) :: ngrid
+!      intent,integer(in) :: nq
+
+!      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
+      real, intent(in) :: pq(ngrid,nlayermx) !water ice mixing ratios (kg/kg)
+      real, intent(in) :: pt(ngrid,nlayermx) !temperature (K)
+      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosol radii
+!      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii
+      real, intent(out) :: nueffrad(ngrid,nlayermx) ! dispersion      
+!      real, intent(out) :: nueffrad(ngrid,nlayermx,naerkind) ! dispersion      
 
       integer :: ig,l
@@ -158,6 +163,8 @@
                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
+!               reffrad(ig,l,iaero_h2o)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
+               reffrad(ig,l)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
+!               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
+               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
             enddo
          enddo
@@ -167,9 +174,13 @@
                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_liq  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o*pi*rhowater) )
+	       zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
+!	       zrad_ice  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+	       zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+!               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
+               nueffrad(ig,l) = 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)
+!	       reffrad(ig,l,iaero_h2o) = min(max(zrad,1.e-6),100.e-6)
+	       reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
                enddo
             enddo      
@@ -193,5 +204,5 @@
 !==================================================================
       use watercommon_h, Only: rhowater,rhowaterice
-      USE tracer_h
+!      USE tracer_h
       Implicit none
 
@@ -201,5 +212,5 @@
 #include "comcstfi.h"
 
-      integer :: ngrid
+      integer,intent(in) :: ngrid
 
       real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
@@ -236,7 +247,7 @@
 !
 !==================================================================
-      use radinc_h, only: naerkind
-      use aerosol_mod, only : iaero_co2
-      USE tracer_h
+!      use radinc_h, only: naerkind
+!      use aerosol_mod, only : iaero_co2
+      USE tracer_h, only:igcm_co2_ice,rho_co2
       Implicit none
 
@@ -246,8 +257,9 @@
 #include "comcstfi.h"
 
-      integer :: ngrid,nq
+      integer,intent(in) :: ngrid,nq
 
       real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
-      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+!      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosols radii (K)
 
       integer :: ig,l
@@ -258,10 +270,12 @@
 
       if (radfixed) then
-         reffrad(1:ngrid,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
+!         reffrad(1:ngrid,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
+         reffrad(1:ngrid,1:nlayermx) = 5.e-5 ! CO2 ice
       else
          do l=1,nlayermx
             do ig=1,ngrid
                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)
+!               reffrad(ig,l,iaero_co2) = min(max(zrad,1.e-6),100.e-6)
+               reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
             enddo
          enddo      
@@ -285,17 +299,19 @@
 !
 !==================================================================
-      use radinc_h, only: naerkind
-      use aerosol_mod, only : iaero_dust
-      Implicit none
-
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-
-      integer :: ngrid
-
-      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+!      use radinc_h, only: naerkind
+!      use aerosol_mod, only : iaero_dust
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+
+!      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosols radii (K)
             
-      reffrad(1:ngrid,1:nlayermx,iaero_dust) = 2.e-6 ! dust
+!      reffrad(1:ngrid,1:nlayermx,iaero_dust) = 2.e-6 ! dust
+      reffrad(1:ngrid,1:nlayermx) = 2.e-6 ! dust
 
    end subroutine dust_reffrad
@@ -315,17 +331,19 @@
 !
 !==================================================================
-      use radinc_h, only: naerkind
-      use aerosol_mod, only : iaero_h2so4
-      Implicit none
-
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-
-      integer :: ngrid
-
-      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+!      use radinc_h, only: naerkind
+!      use aerosol_mod, only : iaero_h2so4
+      Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+
+!      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
+      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosols radii (K)
                 
-      reffrad(1:ngrid,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4
+!      reffrad(1:ngrid,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4
+      reffrad(1:ngrid,1:nlayermx) = 1.e-6 ! h2so4
 
    end subroutine h2so4_reffrad
Index: trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 858)
@@ -28,31 +28,28 @@
 #include "callkeys.h"
 
-      integer ngrid,nq
-
-!     Pre-arguments (for universal model)
-      real pq(ngrid,nlayermx,nq)       ! tracer (kg/kg)
-      real qsurf(ngrid,nq)             ! tracer at the surface (kg.m-2)
-      REAL pdt(ngrid,nlayermx),pdq(ngrid,nlayermx,nq)
-
-      real dqrain(ngrid,nlayermx,nq)   ! tendency of H2O precipitation (kg/kg.s-1)
-      real dqsrain(ngrid)                ! rain flux at the surface (kg.m-2.s-1)
-      real dqssnow(ngrid)                ! snow flux at the surface (kg.m-2.s-1)
-      REAL d_t(ngrid,nlayermx)           ! temperature increment
-
 !     Arguments
-      REAL ptimestep                    ! time interval
-      REAL pplev(ngrid,nlayermx+1)    ! inter-layer pressure
-      REAL pplay(ngrid,nlayermx)      ! mid-layer pressure
-      REAL t(ngrid,nlayermx)          ! input temperature (K)
+      integer,intent(in) :: ngrid ! number of atmospherci columns
+      integer,intent(in) :: nq ! number of tracers
+      real,intent(in) :: ptimestep    ! time interval
+      real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
+      real,intent(in) :: pplay(ngrid,nlayermx)   ! mid-layer pressure (Pa)
+      real,intent(in) :: t(ngrid,nlayermx) ! input temperature (K)
+      real,intent(in) :: pdt(ngrid,nlayermx) ! input tendency on temperature (K/s)      
+      real,intent(in) :: pq(ngrid,nlayermx,nq)  ! tracers (kg/kg)
+      real,intent(in) :: pdq(ngrid,nlayermx,nq) ! input tendency on tracers
+      real,intent(out) :: d_t(ngrid,nlayermx) ! temperature tendency (K/s)
+      real,intent(out) :: dqrain(ngrid,nlayermx,nq) ! tendency of H2O precipitation (kg/kg.s-1)
+      real,intent(out) :: dqsrain(ngrid)  ! rain flux at the surface (kg.m-2.s-1)
+      real,intent(out) :: dqssnow(ngrid)  ! snow flux at the surface (kg.m-2.s-1)
+      real,intent(in) :: rneb(ngrid,nlayermx) ! cloud fraction
+
       REAL zt(ngrid,nlayermx)         ! working temperature (K)
       REAL ql(ngrid,nlayermx)         ! liquid water (Kg/Kg)
       REAL q(ngrid,nlayermx)          ! specific humidity (Kg/Kg)
-      REAL rneb(ngrid,nlayermx)       ! cloud fraction
       REAL d_q(ngrid,nlayermx)        ! water vapor increment
       REAL d_ql(ngrid,nlayermx)       ! liquid water / ice increment
 
 !     Subroutine options
-      REAL seuil_neb                    ! Nebulosity threshold
-      PARAMETER (seuil_neb=0.001)
+      REAL,PARAMETER :: seuil_neb=0.001  ! Nebulosity threshold
 
       INTEGER,save :: precip_scheme      ! id number for precipitaion scheme
@@ -67,13 +64,10 @@
       REAL,SAVE :: c1
 
-      INTEGER ninter
-      PARAMETER (ninter=5)
-
-      logical evap_prec                 ! Does the rain evaporate?
-      parameter(evap_prec=.true.)
+      INTEGER,PARAMETER :: ninter=5
+
+      logical,parameter :: evap_prec=.true. ! Does the rain evaporate?
 
 !     for simple scheme
-      real t_crit
-      PARAMETER (t_crit=218.0)
+      real,parameter :: t_crit=218.0
       real lconvert
 
@@ -100,6 +94,5 @@
       INTEGER, SAVE :: i_ice=0  ! water ice
 
-      LOGICAL firstcall
-      SAVE firstcall
+      LOGICAL,SAVE :: firstcall=.true.
 
 !     Online functions
@@ -108,5 +101,4 @@
       fall2v (zzz) =10.6 * ((zzz)**0.31)  !for use with radii
 
-      DATA firstcall /.true./
 
       IF (firstcall) THEN
Index: trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 853)
+++ trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 858)
@@ -4,5 +4,5 @@
       use comdiurn_h
       USE comgeomfi_h
-      USE tracer_h
+      USE tracer_h, only: igcm_h2o_ice
       implicit none
 
@@ -24,12 +24,12 @@
 #include "callkeys.h"
 
-      integer ngrid, nq
+      integer,intent(in) :: ngrid        ! number of atmospheric columns
+      integer,intent(in) :: nq           ! number of tracers
+      real,intent(in) :: rneb(ngrid,nlayermx)    ! cloud fraction     
+      real,intent(out) :: totalrneb(ngrid)       ! total cloud fraction 
+      real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
+      real,intent(in) :: pq(ngrid,nlayermx,nq)   ! tracers (.../kg_of_air)
+      real,intent(in) :: tau(ngrid,nlayermx)
 
-      real rneb(ngrid,nlayermx)             ! cloud fraction     
-      real pplev(ngrid,nlayermx+1)
-      real pq(ngrid,nlayermx,nq)
-      real tau(ngrid,nlayermx)
-
-      real totalrneb(ngrid)                 ! total cloud fraction 
       real, dimension(nlayermx+1) :: masse
       integer, parameter          :: recovery=7
