Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/aer_sedimnt.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/aer_sedimnt.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/aer_sedimnt.F90	(revision 5202)
@@ -17,5 +17,6 @@
 !-----------------------------------------------------------------------
 
-  USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, f_r_wet, vsed_aer
+  USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, DENSO4B, f_r_wet, f_r_wetB, vsed_aer
+  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
   USE dimphy, ONLY : klon,klev
   USE infotrac_phy
@@ -89,7 +90,13 @@
 
       ! stokes-velocity with cunnigham slip- flow correction
-      ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* &
-         (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK))))
-
+      IF(flag_new_strat_compo) THEN
+         ! stokes-velocity with cunnigham slip- flow correction
+         ZVAER(JL,JK,nb) = 2./9.*(DENSO4B(JL,JK,nb)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wetB(JL,JK,nb)*mdw(nb)/2.)**2.* &
+              (1.+ 2.*zlair(JL,JK)/(f_r_wetB(JL,JK,nb)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wetB(JL,JK,nb)*mdw(nb)/zlair(JL,JK))))
+      ELSE
+         ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* &
+              (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK))))
+      ENDIF
+      
       ZSEDFLX(JL,nb)=ZVAER(JL,JK,nb)*ZRHO
       ZSOLAERB(nb)=ZSOLAERB(nb)+ZDTGDP*ZSEDFLX(JL,nb)
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/aerophys.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/aerophys.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/aerophys.F90	(revision 5202)
@@ -5,23 +5,28 @@
   IMPLICIT NONE
 !
-  REAL,PARAMETER                         :: ropx=1500.0              ! default aerosol particle mass density [kg/m3]
-  REAL,PARAMETER                         :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3]
-  REAL,PARAMETER                         :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3]
-  REAL,PARAMETER                         :: mdwmin=0.002e-6          ! dry diameter of smallest aerosol particles [m]
-  REAL,PARAMETER                         :: V_rat=2.0                ! volume ratio of neighboring size bins
-  REAL,PARAMETER                         :: mfrac_H2SO4=0.75         ! default mass fraction of H2SO4 in the aerosol
-  REAL, PARAMETER                        :: mAIRmol=28.949*1.66E-27  ! Average mass of an air molecule [kg]
-  REAL, PARAMETER                        :: mH2Omol=18.016*1.66E-27  ! Mass of an H2O molecule [kg]
-  REAL, PARAMETER                        :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg]
-  REAL, PARAMETER                        :: mSO2mol=64.06*1.66E-27   ! Mass of an SO2 molecule [kg]
-  REAL, PARAMETER                        :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
-  REAL, PARAMETER                        :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
-  REAL, PARAMETER                        :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
-  REAL, PARAMETER                        :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
-  REAL, PARAMETER                        :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
-  REAL, PARAMETER                        :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
-  REAL, PARAMETER                        :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
-  REAL, PARAMETER                        :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
-  REAL, PARAMETER                        :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg] 
+  REAL,PARAMETER    :: ropx=1500.0              ! default aerosol particle mass density [kg/m3]
+  REAL,PARAMETER    :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3]
+  REAL,PARAMETER    :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3]
+  REAL,PARAMETER    :: mdwmin=0.002e-6          ! dry diameter of smallest aerosol particles [m]
+  REAL,PARAMETER    :: V_rat=2.0                ! volume ratio of neighboring size bins
+  REAL,PARAMETER    :: mfrac_H2SO4=0.75         ! default mass fraction of H2SO4 in the aerosol
+  REAL, PARAMETER   :: mAIRmol=28.949*1.66E-27  ! Average mass of an air molecule [kg]
+  REAL, PARAMETER   :: mH2Omol=18.016*1.66E-27  ! Mass of an H2O molecule [kg]
+  REAL, PARAMETER   :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg]
+  REAL, PARAMETER   :: mSO2mol=64.06*1.66E-27   ! Mass of an SO2 molecule [kg]
+  REAL, PARAMETER   :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
+  REAL, PARAMETER   :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
+  REAL, PARAMETER   :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
+  REAL, PARAMETER   :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
+  REAL, PARAMETER   :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
+  REAL, PARAMETER   :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
+  REAL, PARAMETER   :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
+  REAL, PARAMETER   :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
+  REAL, PARAMETER   :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg]
+  REAL, PARAMETER   :: rgas=8.3145 ! molar gas cste (J⋅K−1⋅mol−1=m3⋅Pa⋅K−1⋅mol−1=kg⋅m2⋅s−2⋅K−1⋅mol−1)
+  ! 
+  REAL, PARAMETER   :: MH2O  =1000.*mH2Omol     ! Mass of 1 molec [g] (18.016*1.66E-24)
+  REAL, PARAMETER   :: MH2SO4=1000.*mH2SO4mol   ! Mass of 1 molec [g] (98.082*1.66E-24)
+  REAL, PARAMETER   :: BOLZ  =1.381E-16         ! Boltzmann constant [dyn.cm/K]
 !
 END MODULE aerophys
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/coagulate.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/coagulate.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/coagulate.F90	(revision 5202)
@@ -26,6 +26,7 @@
   USE aerophys
   USE infotrac_phy
-  USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
-
+  USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB
+  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
+  
   IMPLICIT NONE
 
@@ -43,9 +44,11 @@
   ! local variables in coagulation routine
   INTEGER                                       :: i,j,k,nb,ilon,ilev
-  REAL, DIMENSION(nbtr_bin)                     :: radius ! aerosol particle radius in each bin [m]
+  REAL, DIMENSION(nbtr_bin)                     :: radiusdry ! dry aerosol particle radius in each bin [m]
+  REAL, DIMENSION(nbtr_bin)                     :: radiuswet ! wet aerosol particle radius in each bin [m]
   REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
   REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
   REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
-  REAL, DIMENSION(nbtr_bin)                     :: V    ! Volume of bins
+  REAL, DIMENSION(nbtr_bin)                     :: Vdry ! Volume dry of bins
+  REAL, DIMENSION(nbtr_bin)                     :: Vwet ! Volume wet of bins
   REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
   REAL                                          :: eta  ! Dynamic viscosity of air
@@ -82,15 +85,20 @@
   include "YOMCST.h"
 
-  DO i=1, nbtr_bin
-   radius(i)=mdw(i)/2.
-   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
-  ENDDO
-
-  DO j=1, nbtr_bin
-  DO i=1, nbtr_bin
-   Vij(i,j)= V(i)+V(j)
-  ENDDO
-  ENDDO
-
+! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k
+! just need to be calculated in model initialization because mdw(:) size is fixed
+! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for  
+! dry aerosols
+  DO i=1, nbtr_bin
+     radiusdry(i)=mdw(i)/2.
+     Vdry(i)=radiusdry(i)**3.  !neglecting factor 4*RPI/3
+     Vwet(i)=0.0
+  ENDDO
+
+  DO j=1, nbtr_bin
+     DO i=1, nbtr_bin
+        Vij(i,j)= Vdry(i)+Vdry(j)
+     ENDDO
+  ENDDO
+  
 !--pre-compute the f(i,j,k) from Jacobson equation 13
   ff=0.0 
@@ -100,19 +108,20 @@
     IF (k.EQ.1) THEN
       ff(i,j,k)= 0.0
-    ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
+    ELSEIF (k.GT.1.AND.Vdry(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.Vdry(k)) THEN
       ff(i,j,k)= 1.-ff(i,j,k-1)
     ELSEIF (k.EQ.nbtr_bin) THEN
-      IF (Vij(i,j).GE.v(k)) THEN
+      IF (Vij(i,j).GE.Vdry(k)) THEN
         ff(i,j,k)= 1.
       ELSE
         ff(i,j,k)= 0.0
       ENDIF
-    ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
-      ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
+    ELSEIF (k.LE.(nbtr_bin-1).AND.Vdry(k).LE.Vij(i,j).AND.Vij(i,j).LT.Vdry(k+1)) THEN
+      ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k))
     ENDIF
   ENDDO
   ENDDO
   ENDDO
-
+! End of just need to be calculated at initialization because mdw(:) size is fixed
+  
   DO ilon=1, klon
   DO ilev=1, klev 
@@ -120,9 +129,18 @@
   IF (is_strato(ilon,ilev)) THEN
   !compute actual wet particle radius & volume for every grid box
-  DO i=1, nbtr_bin
-   radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
-   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
-  ENDDO
-
+  IF(flag_new_strat_compo) THEN
+     DO i=1, nbtr_bin
+        radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
+        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
+!!      Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3)
+     ENDDO
+  ELSE
+     DO i=1, nbtr_bin
+        radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
+        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
+!!      Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3)
+     ENDDO
+  ENDIF
+  
 !--Calculations for the coagulation kernel---------------------------------------------------------
 
@@ -150,14 +168,21 @@
   Di=0.0
   DO i=1, nbtr_bin
-   Kn(i)=mnfrpth/radius(i)
-   Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
+      Kn(i)=mnfrpth/radiuswet(i)
+      Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
   ENDDO
 
 !--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
   thvelpar=0.0
-  DO i=1, nbtr_bin
-   m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000.
-   thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
-  ENDDO
+  IF(flag_new_strat_compo) THEN
+     DO i=1, nbtr_bin
+        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000.
+        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
+     ENDDO
+  ELSE
+     DO i=1, nbtr_bin
+        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000.
+        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
+     ENDDO
+  ENDIF
 
 !--pre-compute the particle mean free path mfppar(i) from equation 22
@@ -171,8 +196,9 @@
   delta=0.0
   DO i=1, nbtr_bin
-   delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ &
-           & (6.*radius(i)*mfppar(i))-2.*radius(i)
-  ENDDO
-
+      delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ &
+           & (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i)
+  ENDDO
+
+!   beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j 
 !--pre-compute the beta(i,j) from equation 17 in Jacobson
   num=0.0
@@ -180,20 +206,20 @@
   DO i=1, nbtr_bin
 !
-   num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
-   denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
-        & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
-   beta(i,j)=num/denom
+     num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j))
+     denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
+          & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j)))
+     beta(i,j)=num/denom
 !
 !--compute enhancement factor due to van der Waals forces
    IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
-     Evdw=1.0
+      Evdw=1.0
    ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
-     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
-     xvdW = LOG(1.+AvdWi)
-     EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
+      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
+      xvdW = LOG(1.+AvdWi)
+      EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
    ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
-     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
-     xvdW = LOG(1.+AvdWi)
-     EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
+      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
+      xvdW = LOG(1.+AvdWi)
+      EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
    ENDIF
 !
@@ -209,5 +235,6 @@
   denom=0.0
   DO j=1, nbtr_bin
-  denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
+     !    fraction of coagulation of k and j that is not giving k
+     denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
   ENDDO
 
@@ -219,13 +246,34 @@
     num=0.0
     DO j=1, k
-    numi=0.0
-    DO i=1, k-1
-    numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+       numi=0.0
+       DO i=1, k-1
+!           
+!           see Jacobson: " In order to conserve volume and volume concentration (which
+!           coagulation physically does) while giving up some accuracy in number concentration" 
+!
+!           Coagulation of i and j giving k
+!           with V(i) and then V(j) because it considers i,j and j,i with the double loop
+!
+!           BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation 
+!           kernel beta(i,j) accounts for wet aerosols -> reply below
+!
+!             numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+            numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+       ENDDO
+       num=num+numi
     ENDDO
-    num=num+numi
-    ENDDO
 
 !--calculate new concentration of other bins
-    tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/(1.+pdtcoag*denom)/V(k)
+!      tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
+    tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) )
+!
+!       In constant composition (no dependency on aerosol size because no kelvin effect)
+!       V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i)
+!       so numi and num are proportional (f_r_wet(ilon,ilev)**3)
+!       and so 
+!        tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
+!                     =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) ) 
+!          with num_dry=...beta(i,j)*Vdry(i)*....
+!       so in old STRATAER (.not.flag_new_strat_compo), it was correct
   ENDIF
 
@@ -234,5 +282,5 @@
 !--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
   DO i=1, nbtr_bin
-   tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
+     tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
   ENDDO
 
@@ -240,4 +288,5 @@
   ENDDO !--end of loop klev
   ENDDO !--end of loop klon
+! *********************************************
 
 END SUBROUTINE COAGULATE
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/cond_evap_tstep_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/cond_evap_tstep_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/cond_evap_tstep_mod.F90	(revision 5202)
@@ -9,6 +9,162 @@
 CONTAINS
 
+      SUBROUTINE condens_evapor_rate_kelvin(R2SO4G,t_seri,pplay,R2SO4, &
+          & DENSO4,f_r_wet,R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
+!
+!     INPUT: 
+!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
+!     t_seri: temperature (K)
+!     pplay: pressure (Pa)
+!     R2SO4: aerosol H2SO4 weight fraction (percent) - flat surface (does not depend on aerosol size)
+!     DENSO4: aerosol density (gr/cm3) 
+!     f_r_wet: factor for converting dry to wet radius 
+!        assuming 'flat surface' composition (does not depend on aerosol size)
+!     variables that depends on aerosol size because of Kelvin effect
+!     R2SO4Gik: number density of gaseous H2SO4 [molecules/cm3] - depends on aerosol size
+!     DENSO4ik: aerosol density (gr/cm3) - depends on aerosol size
+!     f_r_wetik: factor for converting dry to wet radius - depends on aerosol size
+!     RRSI: radius [cm]
+
+      USE aerophys
+      USE infotrac_phy
+      USE YOMCST, ONLY : RPI
+      USE sulfate_aer_mod, ONLY : wph2so4, surftension, solh2so4, rpmvh2so4
+      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
+      
+      IMPLICIT NONE
+      
+      REAL, PARAMETER :: third=1./3.
+      
+      ! input variables
+      REAL            :: R2SO4G !H2SO4 number density [molecules/cm3]
+      REAL            :: t_seri
+      REAL            :: pplay
+      REAL            :: R2SO4
+      REAL            :: DENSO4
+      REAL            :: f_r_wet
+      REAL            :: R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
+      
+      ! output variables
+      REAL            :: FL(nbtr_bin)
+      REAL            :: ASO4(nbtr_bin)
+      REAL            :: DNDR(nbtr_bin)
+      
+      ! local variables
+      INTEGER            :: IK
+      REAL            :: ALPHA,CST
+      REAL            :: WH2(nbtr_bin)
+      REAL            :: RP,VTK,AA,FL1,RKNUD
+      REAL            :: DND
+      REAL            :: ATOT,AH2O
+      REAL            :: RRSI_wet(nbtr_bin)
+      REAL            :: FPATH, WPP, XA, FKELVIN
+      REAL            :: surtens, mvh2so4, temp
+      
+! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
+!    ------------------------------------------------------------------
+!                                  EXCEPT CN
+!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
+!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
+
+!                   MOLECULAR ACCOMODATION OF H2SO4
+!     H2SO4 accommodation  coefficient [condensation/evaporation]
+      ALPHA = ALPH2SO4
+!      FPLAIR=(2.281238E-5)*TAIR/PAIR
+!     1.E2 (m to cm), 
+      CST=1.E2*2.281238E-5
+!     same expression as in coagulate
+!     in coagulate: mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
+!     mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
+!     mnfrpth=2.28E-5*t_seri/pplay
+      
+      temp = min( max(t_seri, 190.), 300.) ! 190K <= temp <= 300K
+      
+      RRSI_wet(:)=RRSI(:)*f_r_wetik(:)
+
+!     Pruppa and Klett
+      FPATH=CST*t_seri/pplay
+    
+!     H2SO4 mass fraction in aerosol
+      WH2(:)=R2SO4ik(:)*1.0E-2
+
+!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
+!                               AYERS ET AL (1980)
+!                                  (MU-MU0)
+!      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
+!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
+!      RP=EXP(RP)*1.01325E6/0.086
+!!      RP=EXP(RP)*1.01325E6
+!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
+
+!      DND=RP*6.02E23/(8.31E7*t_seri)
+
+!                                 KELVIN EFFECT FACTOR
+!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
+!!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
+!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)
+
+!                                  MEAN KINETIC VELOCITY
+!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
+!                                  IN CM/SEC
+      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
+!                                 KELVIN EFFECT FACTOR
+
+!     Loop on bin radius (RRSI in cm)
+      DO IK=1,nbtr_bin
+
+      IF(R2SO4ik(IK) > 0.0) THEN
+
+!       h2so4 mass fraction (0<wpp<1)
+        wpp=R2SO4ik(IK)*1.e-2    
+        xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 
+!       equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
+        DND=solh2so4(t_seri,xa)
+!          KELVIN EFFECT:  
+!       surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
+        surtens=surftension(temp,xa)
+!       partial molar volume of h2so4 (cm3.mol-1 =1.e-6.m3.mol-1)
+        mvh2so4= rpmvh2so4(temp,R2SO4ik(IK))
+!       Kelvin factor (MKS)
+        fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2so4/ (1.e-2*RRSI_wet(IK)*rgas*temp) )
+!                              
+        DNDR(IK) =DND*fkelvin
+
+        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))
+
+!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
+!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
+!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
+!       KNUDSEN NUMBER FPATH/RRSI
+!       NEW VERSION (SEE NOTES)
+        RKNUD=FPATH/RRSI_wet(IK)
+!       SENFELD
+        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) & 
+     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
+!       TURCO
+!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
+!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
+!        FL=FL1*RRSI(IK)*RRSI(IK)
+!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )
+
+!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
+        ATOT=4.0*RPI*DENSO4ik(IK)*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
+        ASO4(IK)=WH2(IK)*ATOT/MH2SO4 !attention: g
+!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
+!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
+!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
+        AH2O=(1.0-WH2(IK))*ATOT/MH2O !attention: g
+
+!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
+!       IT IS FOR KEM BUT THERE ARE OTHER WAYS
+
+      ENDIF  
+
+      ENDDO !loop over bins
+
+      END SUBROUTINE condens_evapor_rate_kelvin
+      
+!********************************************************************
       SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, &
-                   & DENSO4,f_r_wet,RRSI,Vbin,FL,ASO4,DNDR)
+                   & DENSO4,f_r_wet,FL,ASO4,DNDR)
 !
 !     INPUT: 
@@ -22,4 +178,5 @@
       USE infotrac_phy
       USE YOMCST, ONLY : RPI
+      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
 
       IMPLICIT NONE
@@ -33,7 +190,5 @@
       REAL DENSO4
       REAL f_r_wet
-      REAL RRSI(nbtr_bin)
-      REAL Vbin(nbtr_bin)
-
+      
       ! output variables
       REAL FL(nbtr_bin)
@@ -48,6 +203,5 @@
       REAL ATOT,AH2O
       REAL RRSI_wet(nbtr_bin)
-      REAL Vbin_wet(nbtr_bin)
-      REAL MH2SO4,MH2O,BOLZ,FPATH
+      REAL FPATH
 
 ! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
@@ -57,20 +211,14 @@
 !     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
 
-!                   WEIGHT OF 1 MOLEC IN G
-      MH2O  =1000.*mH2Omol !18.016*1.66E-24
-      MH2SO4=1000.*mH2SO4mol !98.082*1.66E-24
-!                   BOLTZMANN CONSTANTE IN DYN.CM/K
-      BOLZ  =1.381E-16
 !                   MOLECULAR ACCOMODATION OF H2SO4
-!     raes and van dingen
-      ALPHA =0.1    
+!     H2SO4 accommodation coefficient [condensation/evaporation]
+      ALPHA = ALPH2SO4
 !      FPLAIR=(2.281238E-5)*TAIR/PAIR
 !     1.E2 (m to cm), 
       CST=1.E2*2.281238E-5
 
-      ! compute local wet particle radius and volume
+      ! compute local wet particle radius [cm]
       RRSI_wet(:)=RRSI(:)*f_r_wet
-      Vbin_wet(:)=Vbin(:)*f_r_wet**3
-
+      
 !     Pruppa and Klett
       FPATH=CST*t_seri/pplay
@@ -138,10 +286,11 @@
 
 !********************************************************************
-      SUBROUTINE cond_evap_part(dt,FL,ASO4,f_r_wet,RRSI,Vbin,tr_seri)
+      SUBROUTINE condens_evapor_part(dt,FL,ASO4,f_r_wet,tr_seri)
 
       USE aerophys
       USE infotrac_phy
       USE YOMCST, ONLY : RPI
-
+      USE strataer_local_var_mod, ONLY : RRSI,Vbin
+      
       IMPLICIT NONE
 
@@ -151,10 +300,8 @@
       REAL ASO4(nbtr_bin)
       REAL f_r_wet
-      REAL RRSI(nbtr_bin)
-      REAL Vbin(nbtr_bin)
-
+      
       ! output variables
       REAL tr_seri(nbtr)
-
+      
       ! local variables
       REAL tr_seri_new(nbtr)
@@ -211,5 +358,5 @@
       tr_seri(:)=tr_seri_new(:)
 
-      END SUBROUTINE cond_evap_part
+      END SUBROUTINE condens_evapor_part
 
 END MODULE cond_evap_tstep_mod
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/micphy_tstep.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/micphy_tstep.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/micphy_tstep.F90	(revision 5202)
@@ -8,5 +8,6 @@
   USE aerophys
   USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat
-  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet
+  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, &
+       f_r_wet, R2SO4B, DENSO4B, f_r_wetB
   USE nucleation_tstep_mod
   USE cond_evap_tstep_mod
@@ -14,5 +15,5 @@
   USE YOMCST, ONLY : RPI, RD, RG
   USE print_control_mod, ONLY: lunout
-  USE strataer_local_var_mod
+  USE strataer_local_var_mod ! contains also RRSI and Vbin
   
   IMPLICIT NONE
@@ -35,10 +36,8 @@
   REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
   REAL                      :: x    ! molefraction of H2SO4 in the critical cluster     
-  REAL Vbin(nbtr_bin)
   REAL a_xm, b_xm, c_xm
   REAL PDT, dt
   REAL H2SO4_init
   REAL ACTSO4(klon,klev)
-  REAL RRSI(nbtr_bin)
   REAL nucl_rate
   REAL cond_evap_rate
@@ -48,9 +47,6 @@
   REAL DNDR(nbtr_bin)
   REAL H2SO4_sat
-
-  DO it=1,nbtr_bin
-    Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
-  ENDDO
-
+  REAL R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
+  
   !coefficients for H2SO4 density parametrization used for nucleation if ntot<4
   a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + &
@@ -61,12 +57,9 @@
        & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
 
-  ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
-  CALL STRAACT(ACTSO4)
-
-  ! compute particle radius in cm RRSI from diameter in m 
-  DO it=1,nbtr_bin
-    RRSI(it)=mdw(it)/2.*100.
-  ENDDO
-
+  IF(.not.flag_new_strat_compo) THEN
+     ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
+     CALL STRAACT(ACTSO4)
+  ENDIF
+  
   DO ilon=1, klon
 !
@@ -104,7 +97,16 @@
       ENDIF
       ! compute cond/evap rate in kg(H2SO4)/kgA/s
-      CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
-             & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
-             & RRSI,Vbin,FL,ASO4,DNDR)
+      IF(flag_new_strat_compo) THEN
+         R2SO4ik(:)   = R2SO4B(ilon,ilev,:)
+         DENSO4ik(:)  = DENSO4B(ilon,ilev,:)
+         f_r_wetik(:) = f_r_wetB(ilon,ilev,:)
+         CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+              & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+              & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
+      ELSE
+         CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+              & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+              & FL,ASO4,DNDR)
+      ENDIF
       ! Compute H2SO4 saturate vapor for big particules
       H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol)
@@ -127,7 +129,7 @@
       tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt)
       ! apply cond to bins
-      CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
+      CALL condens_evapor_part(dt,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
       ! apply nucl. to bins
-      CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:))
+      CALL nucleation_part(nucl_rate,ntot,x,dt,tr_seri(ilon,ilev,:))
       ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
       budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
@@ -142,7 +144,13 @@
         & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
     ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys)
-    CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
-           & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
-           & RRSI,Vbin,FL,ASO4,DNDR)
+    IF(flag_new_strat_compo) THEN
+       CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+            & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+            & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
+    ELSE
+       CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+            & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+            & FL,ASO4,DNDR)
+    ENDIF
     ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
     DO it=1,nbtr_bin
@@ -159,5 +167,5 @@
     tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys)
     ! apply evap to bins
-    CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
+    CALL condens_evapor_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
     ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
     budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/miecalc_aer.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/miecalc_aer.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/miecalc_aer.F90	(revision 5202)
@@ -16,5 +16,5 @@
 
   USE phys_local_var_mod, ONLY: tr_seri, mdw, alpha_bin, piz_bin, cg_bin
-  USE aerophys
+  USE aerophys, ONLY: dens_aer_dry, dens_aer_ref, V_rat
   USE aero_mod
   USE infotrac_phy, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat
@@ -226,19 +226,5 @@
     40000.000,    0.2500,   1.48400,   1.0000E-08, &
     50000.000,    0.2000,   1.49800,   1.0000E-08 /), (/nb_lambda_h2so4,4/), order=(/2,1/) )
-
-  !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
-    mdw(1)=mdwmin
-    IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
-      mdw(2)=mdw(1)*2.**(1./3.)
-      DO it=3, nbtr_bin
-        mdw(it)=mdw(it-1)*V_rat**(1./3.)
-      ENDDO 
-    ELSE
-      DO it=2, nbtr_bin
-        mdw(it)=mdw(it-1)*V_rat**(1./3.)
-      ENDDO
-    ENDIF
-    WRITE(lunout,*) 'init mdw=', mdw
-
+     
     !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
     DO bin_number=1, nbtr_bin
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/nucleation_tstep_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/nucleation_tstep_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/nucleation_tstep_mod.F90	(revision 5202)
@@ -70,9 +70,10 @@
 !--------------------------------------------------------------------------------------------------
 
-SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri)
+SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,tr_seri)
 
   USE aerophys
   USE infotrac_phy
-
+  USE strataer_local_var_mod, ONLY : Vbin
+  
   IMPLICIT NONE
 
@@ -82,6 +83,5 @@
   REAL x    ! mole raction of H2SO4 in the critical cluster
   REAL dt
-  REAL Vbin(nbtr_bin)
-
+  
   ! output variable
   REAL tr_seri(nbtr)
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_local_var_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 5202)
@@ -51,4 +51,8 @@
   
   !============= NUCLEATION VARS =============
+  ! MOLECULAR ACCOMODATION OF H2SO4 (Raes and Van Dingen)
+  REAL,SAVE    :: ALPH2SO4               ! H2SO4 accommodation  coefficient [condensation/evaporation]
+  !$OMP THREADPRIVATE(ALPH2SO4)
+  
   ! flag to constraint nucleation rate in a lat/pres box
   LOGICAL,SAVE :: flag_nuc_rate_box      ! Nucleation rate limit or not to a lat/pres
@@ -64,6 +68,4 @@
   INTEGER,SAVE :: flh2o  ! ds stratemit : flh2o =0 (tr_seri), flh2o=1 (dq) 
   !$OMP THREADPRIVATE(flh2o)
-!  REAL,ALLOCATABLE,SAVE    :: d_q_emiss(:,:)
-!  !$OMP THREADPRIVATE(d_q_emiss)
   
   REAL,ALLOCATABLE,SAVE    :: budg_emi(:,:)            !DIMENSION(klon,n)
@@ -144,4 +146,7 @@
   !$OMP THREADPRIVATE(day_emit_roc)
   
+  REAL,ALLOCATABLE,SAVE    :: RRSI(:) ! radius [cm] for each aerosol size
+  REAL,ALLOCATABLE,SAVE    :: Vbin(:) ! volume [m3] for each aerosol size  
+  !$OMP THREADPRIVATE(RRSI, Vbin)
   REAL,SAVE    :: dlat, dlon             ! delta latitude and d longitude of grid in degree
   !$OMP THREADPRIVATE(dlat, dlon)
@@ -153,5 +158,10 @@
     USE print_control_mod, ONLY : lunout
     USE mod_phys_lmdz_para, ONLY : is_master
-    USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas
+    USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas,nbtr_bin
+    USE phys_local_var_mod, ONLY : mdw
+    USE aerophys, ONLY: mdwmin, V_rat
+    USE YOMCST  , ONLY : RPI
+    
+    INTEGER :: it
     
     WRITE(lunout,*) 'IN STRATAER_LOCAL_VAR INIT WELCOME!'
@@ -185,4 +195,5 @@
     
     ! nuc init
+    ALPH2SO4 = 0.1
     flag_nuc_rate_box = .FALSE.
     nuclat_min=0  ; nuclat_max=0
@@ -238,4 +249,33 @@
     ENDIF ! if master
     
+    !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
+    mdw(1)=mdwmin
+    IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
+       mdw(2)=mdw(1)*2.**(1./3.)
+       DO it=3, nbtr_bin
+          mdw(it)=mdw(it-1)*V_rat**(1./3.)
+       ENDDO
+    ELSE
+       DO it=2, nbtr_bin
+          mdw(it)=mdw(it-1)*V_rat**(1./3.)
+       ENDDO
+    ENDIF
+    IF (is_master) WRITE(lunout,*) 'init mdw=', mdw
+    
+    !   compute particle radius RRSI [cm] and volume Vbin [m3] from diameter mdw [m]
+    ALLOCATE(RRSI(nbtr_bin), Vbin(nbtr_bin))
+    
+    DO it=1,nbtr_bin
+       !     [cm]
+       RRSI(it)=mdw(it)/2.*100.
+       !     [m3]
+       Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
+    ENDDO
+    
+    IF (is_master) THEN
+       WRITE(lunout,*) 'init RRSI=', RRSI
+       WRITE(lunout,*) 'init Vbin=', Vbin
+    ENDIF
+    
     WRITE(lunout,*) 'IN STRATAER INIT END'
     
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_nuc_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_nuc_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/strataer_nuc_mod.F90	(revision 5202)
@@ -13,5 +13,6 @@
     USE print_control_mod, ONLY : lunout
     USE mod_phys_lmdz_para, ONLY : is_master
-    USE strataer_local_var_mod, ONLY: flag_nuc_rate_box,nuclat_min,nuclat_max,nucpres_min,nucpres_max
+    USE strataer_local_var_mod, ONLY: ALPH2SO4,flag_nuc_rate_box,nuclat_min,nuclat_max, &
+         nucpres_min,nucpres_max
     
     !Config Key  = flag_nuc_rate_box
@@ -30,6 +31,10 @@
     CALL getin_p('nucpres_max',nucpres_max)
     
+    ! Read argument H2SO4 accommodation  coefficient [condensation/evaporation]
+    CALL getin_p('alph2so4',ALPH2SO4)
+    
     !============= Print params =============
     IF (is_master) THEN
+       WRITE(lunout,*) 'IN STRATAER_NUC : ALPH2SO4 = ',alph2so4
        WRITE(lunout,*) 'IN STRATAER_NUC : flag_nuc_rate_box = ',flag_nuc_rate_box
        IF (flag_nuc_rate_box) THEN
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/sulfate_aer_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/sulfate_aer_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/sulfate_aer_mod.F90	(revision 5202)
@@ -7,73 +7,145 @@
 
 !*******************************************************************
-  SUBROUTINE STRACOMP_BIN(sh,t_seri,pplay)
-!
-!   Aerosol H2SO4 weight fraction as a function of PH2O and temperature
-!   INPUT:
-!   sh: VMR of H2O
-!   t_seri: temperature (K)
-!   pplay: middle layer pression (Pa)
-!
-!   OUTPUT:
-!   R2SO4: aerosol H2SO4 weight fraction (percent)
+      SUBROUTINE STRACOMP_KELVIN(sh,t_seri,pplay)
+!
+!     Aerosol H2SO4 weight fraction as a function of PH2O and temperature
+!     INPUT:
+!     sh: MMR of H2O
+!     t_seri: temperature (K)
+!     pplay: middle layer pression (Pa)
+!
+!     Modified in modules:
+!     R2SO4: aerosol H2SO4 weight fraction (percent)
+!     R2SO4B: aerosol H2SO4 weight fraction (percent) for each aerosol bin
+!     DENSO4: aerosol density (gr/cm3)
+!     DENSO4B: aerosol density (gr/cm3)for each aerosol bin
+!     f_r_wet: factor for converting dry to wet radius 
+!        assuming 'flat surface' composition (does not depend on aerosol size)
+!     f_r_wetB: factor for converting dry to wet radius 
+!        assuming 'curved surface' composition (depends on aerosol size)
     
-    USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands
-    USE aerophys
-    USE phys_local_var_mod, ONLY: R2SO4
+      USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands
+      USE infotrac_phy, ONLY : nbtr_bin
+      USE aerophys
+      USE phys_local_var_mod, ONLY: R2SO4, R2SO4B, DENSO4, DENSO4B, f_r_wet, f_r_wetB
+      USE strataer_local_var_mod, ONLY: RRSI
+!     WARNING: in phys_local_var_mod R2SO4B, DENSO4B, f_r_wetB (klon,klev,nbtr_bin)
+!          and dens_aer_dry must be declared somewhere
     
-    IMPLICIT NONE
+      IMPLICIT NONE
     
-    REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
-    REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression in the middle of each layer (Pa)
-    REAL,DIMENSION(klon,klev),INTENT(IN)          :: sh      ! specific humidity
+      REAL,DIMENSION(klon,klev),INTENT(IN)    :: t_seri  ! Temperature
+      REAL,DIMENSION(klon,klev),INTENT(IN)    :: pplay   ! pression in the middle of each layer (Pa)
+      REAL,DIMENSION(klon,klev),INTENT(IN)    :: sh      ! specific humidity (kg h2o/kg air)
+      
+!     local variables
+      integer         :: ilon,ilev,ik
+      real, parameter :: rath2oair = mAIRmol/mH2Omol  
+      real, parameter :: third = 1./3.
+      real            :: pph2ogas(klon,klev) 
+      real            :: temp, wpp, xa, surtens, mvh2o, radwet, fkelvin, pph2okel, r2so4ik, denso4ik
+!----------------------------------------
+  
+!     gas-phase h2o partial pressure (Pa) 
+!                                vmr=sh*rath2oair
+      pph2ogas(:,:) = pplay(:,:)*sh(:,:)*rath2oair
     
-    REAL ks(7)
-    REAL t,qh2o,ptot,pw
-    REAL a,b,c,det
-    REAL xsb,msb
+      DO ilon=1,klon
+      DO ilev=1,klev
+          
+        temp = max(t_seri(ilon,ilev),190.)
+        temp = min(temp,300.)
+
+!    ***   H2SO4-H2O flat surface ***
+!!       equilibrium H2O pressure over pure flat liquid water (Pa)
+!!        pflath2o=psh2o(temp)
+!       h2so4 weight percent(%) = f(P_h2o(Pa),T)
+        R2SO4(ilon,ilev)=wph2so4(pph2ogas(ilon,ilev),temp)  
+!       h2so4 mass fraction (0<wpp<1)
+        wpp=R2SO4(ilon,ilev)*1.e-2     
+!       mole fraction
+        xa=18.*wpp/(18.*wpp+98.*(1.-wpp))
+
+!        CHECK:compare h2so4 sat/ pressure (see Marti et al., 97 & reef. therein)
+!               R2SO4(ilon,ilev)=70.    temp=298.15
+!        equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
+!        include conversion from molec/cm3 to Pa
+!        ph2so4=solh2so4(temp,xa)*(1.38065e-16*temp)/10.
+!        print*,' ph2so4=',ph2so4,temp,R2SO4(ilon,ilev)
+!        good match with Martin, et Ayers, not with Gmitro (the famous 0.086)
+
+!       surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
+        surtens=surftension(temp,xa)
+!       molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1)
+        mvh2o= rmvh2o(temp)
+!       aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
+        DENSO4(ilon,ilev)=density(temp,wpp) 
+!           ->x1000., to have it in kg/m3
+!       factor for converting dry to wet radius
+        f_r_wet(ilon,ilev) = (dens_aer_dry/(DENSO4(ilon,ilev)*1.e3)/ & 
+                   &    (R2SO4(ilon,ilev)*1.e-2))**third
+!    ***   End of H2SO4-H2O flat surface ***
+
+
+!          Loop on bin radius (RRSI in cm)
+           DO IK=1,nbtr_bin 
+ 
+!      ***   H2SO4-H2O curved surface - Kelvin effect factor ***
+!            wet radius (m) (RRSI(IK) in [cm])
+             if (f_r_wetB(ilon,ilev,IK) .gt. 1.0) then 
+               radwet = 1.e-2*RRSI(IK)*f_r_wetB(ilon,ilev,IK)
+             else 
+!              H2SO4-H2O flat surface, only on the first timestep
+               radwet = 1.e-2*RRSI(IK)*f_r_wet(ilon,ilev)
+             endif
+!            Kelvin factor: 
+!            surface tension (mN/m=1.e-3.kg/s2) 
+!            molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1)
+             fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2o/ (radwet*rgas*temp) )
+!            equilibrium: pph2o(gas) = pph2o(liq) = pph2o(liq_flat) * fkelvin
+!            equilibrium: pph2o(liq_flat) = pph2o(gas) / fkelvin
+!            h2o liquid partial pressure before Kelvin effect (Pa)
+             pph2okel = pph2ogas(ilon,ilev) / fkelvin
+!            h2so4 weight percent(%) = f(P_h2o(Pa),temp)
+             r2so4ik=wph2so4(pph2okel,temp)
+!            h2so4 mass fraction (0<wpp<1)
+             wpp=r2so4ik*1.e-2    
+!            mole fraction
+             xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 
+!            aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
+             denso4ik=density(temp,wpp)
+!           
+!            recalculate Kelvin factor with surface tension and radwet 
+!                              with new R2SO4B and DENSO4B
+             surtens=surftension(temp,xa)
+!            wet radius (m)
+             radwet = 1.e-2*RRSI(IK)*(dens_aer_dry/(denso4ik*1.e3)/ & 
+                   &    (r2so4ik*1.e-2))**third
+             fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2o / (radwet*rgas*temp) ) 
+             pph2okel=pph2ogas(ilon,ilev) / fkelvin
+!            h2so4 weight percent(%) = f(P_h2o(Pa),temp)
+             R2SO4B(ilon,ilev,IK)=wph2so4(pph2okel,temp)
+!            h2so4 mass fraction (0<wpp<1)
+             wpp=R2SO4B(ilon,ilev,IK)*1.e-2    
+             xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 
+!            aerosol density (gr/cm3) = f(T,h2so4 mass fraction)
+             DENSO4B(ilon,ilev,IK)=density(temp,wpp)
+!            factor for converting dry to wet radius
+             f_r_wetB(ilon,ilev,IK) = (dens_aer_dry/(DENSO4B(ilon,ilev,IK)*1.e3)/ & 
+                   &    (R2SO4B(ilon,ilev,IK)*1.e-2))**third
+!
+!             print*,'R,Rwet(m),kelvin,h2so4(%),ro=',RRSI(ik),radwet,fkelvin, &
+!              &  R2SO4B(ilon,ilev,IK),DENSO4B(ilon,ilev,IK)
+!             print*,' equil.h2so4(molec/cm3), & 
+!              & sigma',solh2so4(temp,xa),surftension(temp,xa)
+
+           ENDDO
+
+      ENDDO
+      ENDDO
+
+      RETURN
     
-    INTEGER ilon,ilev
-    DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/
-    
-!*******************************************************************
-!***     liquid aerosols process
-!*******************************************************************
-!        BINARIES LIQUID AEROROLS:
-    
-    DO ilon=1,klon
-       DO ilev=1,klev
-          
-          t = max(t_seri(ilon,ilev),185.)
-          qh2o=sh(ilon,ilev)/18.*28.9
-          ptot=pplay(ilon,ilev)/100.
-          pw = qh2o*ptot/1013.0
-          pw = min(pw,2.e-3/1013.)
-          pw = max(pw,2.e-5/1013.)
-          
-!*******************************************************************
-!***     binaries aerosols h2so4/h2o
-!*******************************************************************
-          a = ks(3) + ks(4)/t
-          b = ks(1) + ks(2)/t
-          c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw)
-          
-          det = b**2 - 4.*a*c
-          
-          IF (det > 0.) THEN
-             xsb = (-b - sqrt(det))/(2.*a)
-             msb = 55.51*xsb/(1.0 - xsb)
-          ELSE
-             msb = 0.
-          ENDIF
-          R2SO4(ilon,ilev) = 100*msb*0.098076/(1.0 + msb*0.098076)
-          
-          ! H2SO4 min dilution: 0.5%
-          R2SO4(ilon,ilev) = max( R2SO4(ilon,ilev), 0.005 )
-       ENDDO
-    ENDDO
-100 RETURN
-    
-  END SUBROUTINE STRACOMP_BIN
-
+  END SUBROUTINE STRACOMP_KELVIN
 !********************************************************************
     SUBROUTINE STRACOMP(sh,t_seri,pplay)
@@ -544,76 +616,4 @@
 
     END SUBROUTINE
-
-!****************************************************************
-    SUBROUTINE DENH2SA_TABA(t_seri)
-
-!   AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
-!   from Tabazadeh et al. (1994) abaques
-!   ---------------------------------------------
-
-!   
-!   INPUT: 
-!   R2SO4: aerosol H2SO4 weight fraction (percent)
-!   t_seri: temperature (K)
-!   klon: number of latitude bands in the model domain
-!   klev: number of altitude bands in the model domain
-!   for IFS: perhaps add another dimension for longitude
-!
-!   OUTPUT: 
-!   DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
-!   
-    USE dimphy, ONLY : klon,klev
-    USE phys_local_var_mod, ONLY: R2SO4, DENSO4
-    
-    IMPLICIT NONE
-    
-    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
-        
-    INTEGER i,j
-    
-!----------------------------------------------------------------------
-!	... Local variables
-!----------------------------------------------------------------------
-      real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3
-      
-      real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8
-      real :: c1, c2, c3, c4, w
-      
-      
-!   Loop on model domain (2 dimension for UPMC model; 3 for IFS)
-    DO i=1,klon
-       DO j=1,klev
-!----------------------------------------------------------------------
-!	... Temperature variables
-!----------------------------------------------------------------------
-          c1 = t_seri(I,J)- 273.15
-          c2 = c1**2
-          c3 = c1*c2
-          c4 = c1*c3
-!----------------------------------------------------------------------
-!	Polynomial Coefficients
-!----------------------------------------------------------------------
-          a0 = 999.8426 + 334.5402e-4*c1 - 569.1304e-5*c2
-          a1 = 547.2659 - 530.0445e-2*c1 + 118.7671e-4*c2 + 599.0008e-6*c3
-          a2 = 526.295e1 + 372.0445e-1*c1 + 120.1909e-3*c2 - 414.8594e-5*c3 + 119.7973e-7*c4
-          a3 = -621.3958e2 - 287.7670*c1 - 406.4638e-3*c2 + 111.9488e-4*c3 + 360.7768e-7*c4
-          a4 = 409.0293e3 + 127.0854e1*c1 + 326.9710e-3*c2 - 137.7435e-4*c3 - 263.3585e-7*c4
-          a5 = -159.6989e4 - 306.2836e1*c1 + 136.6499e-3*c2 + 637.3031e-5*c3
-          a6 = 385.7411e4 + 408.3717e1*c1 - 192.7785e-3*c2
-          a7 = -580.8064e4 - 284.4401e1*c1
-          a8 = 530.1976e4 + 809.1053*c1
-!----------------------------------------------------------------------
-!	... Summation
-!----------------------------------------------------------------------
-!     w : H2SO4 Weight fraction
-          w=r2SO4(i,j)*0.01
-          DENSO4(i,j) = 0.001*(a0 + w*(a1 + w*(a2 + w*(a3 + w*(a4 +  & 
-               w*(a5 + w*(a6 + w*(a7 + w*(a8 + w*(a9 + w*a10))))))))))
-          DENSO4(i,j) = max (0.0, DENSO4(i,j) )
-
-       ENDDO
-    ENDDO
-
-  END SUBROUTINE DENH2SA_TABA
   
 !****************************************************************
@@ -764,4 +764,299 @@
        RETURN
        END SUBROUTINE
-
+!********************************************************************
+!-----------------------------------------------------------------------
+      real function psh2so4(T) result(psh2so4_out)
+!     equilibrium H2SO4 pressure over pure H2SO4 solution (Pa)
+!
+!---->Ayers et.al. (1980), GRL (7) pp 433-436
+!     plus corrections for lower temperatures by Kulmala and Laaksonen (1990)
+!     and Noppel et al. (1990)
+
+      implicit none
+      real, intent(in) :: T
+      real, parameter ::      &
+              &  b1=1.01325e5, &
+              &  b2=11.5,  &
+              &  b3=1.0156e4,  &
+              &  b4=0.38/545., &
+              &  tref=360.15
+
+!     saturation vapor pressure ( N/m2 = Pa = kg/(m.s2) )
+      psh2so4_out=b1*exp(  -b2 +b3*( 1./tref-1./T  &
+           &  +b4*(1.+log(tref/T)-tref/T) )   )  
+
+       return
+    end function psh2so4
+!-----------------------------------------------------------------------
+    real function ndsh2so4(T) result(ndsh2so4_out)
+!     equilibrium H2SO4 number density over pure H2SO4 (molec/cm3)
+
+      implicit none
+      real, intent(in) :: T
+      real :: presat
+
+!     Boltzmann constant ( 1.38065e-23 J/K = m2⋅kg/(s2⋅K) )
+!      akb idem in cm2⋅g/(s2⋅K) 
+      real, parameter :: akb=1.38065e-16
+
+!     pure h2so4 saturation vapor pressure (Pa)
+      presat=psh2so4(T)
+!     saturation number density (1/cm3) - (molec/cm3)
+      ndsh2so4_out=presat*10./(akb*T)
+
+       return
+     end function ndsh2so4
+!-----------------------------------------------------------------------
+     real function psh2o(T) result(psh2o_out)
+!     equilibrium H2O pressure over pure liquid water (Pa)
+!
+      implicit none
+      real, intent(in) :: T
+
+      if(T.gt.229.) then
+!        Preining et al., 1981 (from Kulmala et al., 1998)
+!        saturation vapor pressure (N/m2 = 1 Pa = 1 kg/(m·s2))
+         psh2o_out=exp( 77.34491296  -7235.424651/T &
+             &                 -8.2*log(T) + 5.7133e-3*T )
+      else
+!        Tabazadeh et al., 1997, parameterization for 185<T<260
+!        saturation water vapor partial pressure (mb = hPa =1.E2 kg/(m·s2))
+!        or from Clegg and Brimblecombe , J. Chem. Eng., p43, 1995.
+;
+         psh2o_out=18.452406985 -3505.1578807/T &
+              &    -330918.55082/(T*T)             &
+              &    +12725068.262/(T*T*T)
+!        in Pa
+         psh2o_out=100.*exp(psh2o_out)
+      end if
+!      print*,psh2o_out
+      
+       return
+     end function psh2o
+!-----------------------------------------------------------------------
+     real function density(T,so4mfrac) result(density_out)
+!        calculation of particle density (gr/cm3)
+
+!        requires Temperature (T) and acid mass fraction (so4mfrac)
+!---->Vehkamaeki et al. (2002)
+
+      implicit none
+      real, intent(in) :: T, so4mfrac
+      real, parameter :: &
+           &      a1= 0.7681724,&
+           &      a2= 2.184714, &
+           &      a3= 7.163002, &
+           &      a4=-44.31447, &
+           &      a5= 88.74606, &
+           &      a6=-75.73729, &
+           &      a7= 23.43228
+      real, parameter :: &
+           &      b1= 1.808225e-3, &
+           &      b2=-9.294656e-3, &
+           &      b3=-3.742148e-2, &
+           &      b4= 2.565321e-1, &
+           &      b5=-5.362872e-1, &
+           &      b6= 4.857736e-1, &
+           &      b7=-1.629592e-1
+      real, parameter :: &
+           &      c1=-3.478524e-6, &
+           &      c2= 1.335867e-5, &
+           &      c3= 5.195706e-5, &
+           &      c4=-3.717636e-4, &
+           &      c5= 7.990811e-4, &
+           &      c6=-7.458060e-4, &
+           &      c7= 2.581390e-4
+      real :: a,b,c,so4m2,so4m3,so4m4,so4m5,so4m6
+      
+      so4m2=so4mfrac*so4mfrac
+      so4m3=so4mfrac*so4m2
+      so4m4=so4mfrac*so4m3
+      so4m5=so4mfrac*so4m4
+      so4m6=so4mfrac*so4m5
+
+      a=+a1+a2*so4mfrac+a3*so4m2+a4*so4m3 &
+         &     +a5*so4m4+a6*so4m5+a7*so4m6
+      b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3 &
+         &     +b5*so4m4+b6*so4m5+b7*so4m6
+      c=+c1+c2*so4mfrac+c3*so4m2+c4*so4m3 &
+         &     +c5*so4m4+c6*so4m5+c7*so4m6
+      density_out=(a+b*T+c*T*T) ! units are gm/cm**3
+
+       return
+     end function density
+!-----------------------------------------------------------------------
+     real function surftension(T,so4frac) result(surftension_out)
+!        calculation of surface tension (mN/meter)
+!        requires Temperature (T) and acid mole fraction (so4frac)
+!---->Vehkamaeki et al. (2002)
+
+      implicit none
+      real,intent(in) :: T, so4frac
+      real :: a,b,so4mfrac,so4m2,so4m3,so4m4,so4m5,so4sig
+      real, parameter :: &
+            &     a1= 0.11864, &
+            &     a2=-0.11651, &
+            &     a3= 0.76852, &
+            &     a4=-2.40909, &
+            &     a5= 2.95434, &
+            &     a6=-1.25852
+      real, parameter :: &
+            &     b1=-1.5709e-4, &
+            &     b2= 4.0102e-4, &
+            &     b3=-2.3995e-3, &
+            &     b4= 7.611235e-3, &
+            &     b5=-9.37386e-3, &
+            &     b6= 3.89722e-3
+      real, parameter :: convfac=1.e3  ! convert from newton/m to dyne/cm
+      real, parameter :: Mw=18.01528, Ma=98.079
+
+!     so4 mass fraction
+      so4mfrac=Ma*so4frac/( Ma*so4frac+Mw*(1.-so4frac) )
+      so4m2=so4mfrac*so4mfrac
+      so4m3=so4mfrac*so4m2
+      so4m4=so4mfrac*so4m3
+      so4m5=so4mfrac*so4m4
+
+      a=+a1+a2*so4mfrac+a3*so4m2+a4*so4m3+a5*so4m4+a6*so4m5
+      b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3+b5*so4m4+b6*so4m5
+      so4sig=a+b*T
+      surftension_out=so4sig*convfac
+
+       return
+     end function surftension
+!-----------------------------------------------------------------------
+     real function wph2so4(pph2o,T) result(wph2so4_out)
+!     Calculates the equilibrium composition of h2so4 aerosols
+!     as a function of temperature and  H2O pressure, using
+!     the parameterization of Tabazadeh et al., GRL, p1931, 1997.
+!
+!   Parameters
+!
+!    input:
+!      T.....temperature (K)
+!      pph2o..... amhbiant 2o pressure (Pa)
+!
+!    output:
+!      wph2so4......sulfuric acid composition (weight percent wt % h2so4)
+!                     = h2so4 mass fraction*100.
+!
+      implicit none
+      real, intent(in) :: pph2o, T 
+      
+      real :: aw, rh, y1, y2, sulfmolal
+  
+!       psh2o(T): equilibrium H2O pressure over pure liquid water (Pa)
+!       relative humidity 
+        rh=pph2o/psh2o(T)
+!       water activity
+!        aw=min( 0.999,max(1.e-3,rh) )
+        aw=min( 0.999999999,max(1.e-8,rh) ) 
+
+!       composition
+!       calculation of h2so4 molality
+            if(aw .le. 0.05 .and. aw .gt. 0.) then
+               y1=12.372089320*aw**(-0.16125516114) &
+                 &  -30.490657554*aw -2.1133114241
+               y2=13.455394705*aw**(-0.19213122550) &
+                 &  -34.285174607*aw -1.7620073078
+            else if(aw .le. 0.85 .and. aw .gt. 0.05) then
+               y1=11.820654354*aw**(-0.20786404244) &
+                 &  -4.8073063730*aw -5.1727540348
+               y2=12.891938068*aw**(-0.23233847708) &
+                 &  -6.4261237757*aw -4.9005471319
+            else
+               y1=-180.06541028*aw**(-0.38601102592) &
+                 &  -93.317846778*aw +273.88132245
+               y2=-176.95814097*aw**(-0.36257048154) &
+                 &  -90.469744201*aw +267.45509988
+            end if
+!        h2so4 molality (m=moles of h2so4 (solute)/ kg of h2o(solvent))
+         sulfmolal = y1+((T-190.)*(y2-y1)/70.)
+
+!        for a solution containing mh2so4 and mh2o: 
+!        sulfmolal = (mh2so4(gr)/h2so4_molar_mass(gr/mole)) / (mh2o(gr)*1.e-3)
+!        mh2o=1.e3*(mh2so4/Mh2so4)/sulfmolal=1.e3*mh2so4/(Mh2so4*sulfmolal)
+!        h2so4_mass_fraction = mfh2so4 = mh2so4/(mh2o + mh2so4)
+!        mh2o=mh2so4*(1-mfh2so4)/mfh2so4
+!        combining the 2 equations
+!        1.e3*mh2so4/(Mh2so4*sulfmolal) = mh2so4*(1-mfh2so4)/mfh2so4
+!        1.e3/(Mh2so4*sulfmolal) = (1-mfh2so4)/mfh2so4
+!        1000*mfh2so4 = (1-mfh2so4)*Mh2so4*sulfmolal
+!        mfh2so4*(1000.+Mh2so4*sulfmolal) = Mh2so4*sulfmolal
+!        mfh2so4 = Mh2so4*sulfmolal / (1000.+Mh2so4*sulfmolal) 
+!        wph2so4 (% mass fraction)= 100.*Mh2so4*sulfmolal / (1000.+Mh2so4*sulfmolal) 
+!        recall activity of i = a_i = P_i/P_pure_i and 
+!          activity coefficient of i = gamma_i = a_i/X_i (X_i: mole fraction of i)
+!        so  P_i = gamma_i*X_i*P_pure_i
+!        if ideal solution, gamma_i=1, P_i = X_i*P_pure_i
+
+!        h2so4 weight precent
+         wph2so4_out = 9800.*sulfmolal/(98.*sulfmolal+1000.)
+!         print*,rh,pph2o,psh2o(T),vpice(T)
+!         print*,T,aw,sulfmolal,wph2so4_out
+         wph2so4_out = max(wph2so4_out,15.)
+         wph2so4_out = min(wph2so4_out,99.999)
+
+       return
+     end function wph2so4
+!-----------------------------------------------------------------------
+     real function solh2so4(T,xa) result(solh2so4_out)
+!     equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
+
+      implicit none
+      real, intent(in) :: T, xa       ! T(K)  xa(H2SO4 mass fraction)
+      
+      real :: xw, a12,b12, cacta, presat
+      
+      xw=1.0-xa
+
+!     pure h2so4 saturation number density (molec/cm3)
+      presat=ndsh2so4(T)
+!     compute activity of acid 
+      a12=5.672E3 -4.074E6/T +4.421E8/(T*T)
+      b12=1./0.527
+      cacta=10.**(a12*xw*xw/(xw+b12*xa)**2/T)
+!     h2so4 saturation number density over H2SO4/H2O solution (molec/cm3)
+      solh2so4_out=cacta*xa*presat
+
+       return
+     end function solh2so4
+!-----------------------------------------------------------------------      
+     real function rpmvh2so4(T,ws) result(rpmvh2so4_out)
+!     partial molar volume of h2so4 in h2so4/h2o solution (cm3/mole)
+
+      implicit none
+      real, intent(in) :: T, ws
+      real, dimension(22),parameter :: x=(/  &
+       & 2.393284E-02,-4.359335E-05,7.961181E-08,0.0,-0.198716351, &
+       & 1.39564574E-03,-2.020633E-06,0.51684706,-3.0539E-03,4.505475E-06, &
+       & -0.30119511,1.840408E-03,-2.7221253742E-06,-0.11331674116, &
+       & 8.47763E-04,-1.22336185E-06,0.3455282,-2.2111E-03,3.503768245E-06, &
+       & -0.2315332,1.60074E-03,-2.5827835E-06/)
+      
+      real :: w
+
+        w=ws*0.01
+        rpmvh2so4_out=x(5)+x(6)*T+x(7)*T*T+(x(8)+x(9)*T+x(10)*T*T)*w &
+          +(x(11)+x(12)*T+x(13)*T*T)*w*w
+!       h2so4 partial molar volume in h2so4/h2o solution (cm3/mole)
+        rpmvh2so4_out=rpmvh2so4_out*1000.
+        
+       return
+     end function rpmvh2so4
+!-----------------------------------------------------------------------
+     real function rmvh2o(T) result(rmvh2o_out)
+!     molar volume of pure h2o (cm3/mole)
+
+       implicit none
+       real, intent(in) :: T 
+       real, parameter :: x1=2.393284E-02,x2=-4.359335E-05,x3=7.961181E-08
+
+!      1000: L/mole ->  cm3/mole
+!      pure h2o molar volume (cm3/mole)
+       rmvh2o_out=(x1+x2*T+x3*T*T)*1000.
+       
+       return
+     end function rmvh2o
+!
 END MODULE sulfate_aer_mod
Index: LMDZ6/branches/cirrus/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/StratAer/traccoag_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/StratAer/traccoag_mod.F90	(revision 5202)
@@ -9,8 +9,9 @@
        presnivs, xlat, xlon, pphis, pphi, &
        t_seri, pplay, paprs, sh, rh, tr_seri)
-
+    
     USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & 
-        & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part 
-
+        & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
+        & R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode, reff_sulfate
+    
     USE dimphy
     USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
@@ -56,5 +57,5 @@
     REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
     REAL                                   :: m_aer               ! aerosol mass
-    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
+    INTEGER                                :: it, k, i, j, ilon, ilev, itime, i_int, ieru
     LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
     REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
@@ -82,4 +83,5 @@
     INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
     INTEGER                                :: yr, is_bissext
+    REAL                                   :: samoment2, samoment3! 2nd and 3rd order moments of size distribution
 
     IF (is_mpi_root .AND. flag_verbose_strataer) THEN
@@ -88,4 +90,5 @@
     ENDIF
     
+    !   radius [m]
     DO it=1, nbtr_bin
       r_bin(it)=mdw(it)/2.
@@ -117,13 +120,9 @@
 
     IF(flag_new_strat_compo) THEN
-       IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Tabazadeh 1994', flag_new_strat_compo
-       ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) : binary routine (from reprobus)
-       ! H2SO4 mass fraction in aerosol (%) from Tabazadeh et al. (1994).
-       CALL stracomp_bin(sh,t_seri,pplay)
-       
-       ! aerosol density (gr/cm3) - from Tabazadeh
-       CALL denh2sa_taba(t_seri)
+       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
+       ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
+       CALL stracomp_kelvin(sh,t_seri,pplay)
     ELSE
-       IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Bekki 2D model', flag_new_strat_compo
+       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
        ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) 
        ! H2SO4 mass fraction in aerosol (%)
@@ -132,9 +131,9 @@
        ! aerosol density (gr/cm3)
        CALL denh2sa(t_seri)
+       
+       ! compute factor for converting dry to wet radius (for every grid box)
+       f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
     ENDIF
     
-! compute factor for converting dry to wet radius (for every grid box)
-    f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
-
 !--calculate mass of air in every grid box
     DO ilon=1, klon
@@ -348,4 +347,68 @@
     ENDDO
     
+!--compute 
+!     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
+!     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
+!     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
+!     nd_mode: DRY(?) particle concentration in different modes (part/m3)
+     sulfmmr(:,:)=0.0
+     SAD_sulfate(:,:)=0.0
+     sulfmmr_mode(:,:,:)=0.0
+     nd_mode(:,:,:)=0.0
+     reff_sulfate(:,:)=0.0
+     
+     DO i=1,klon
+        DO j=1,klev
+           samoment2=0.0
+           samoment3=0.0
+           DO it=1, nbtr_bin
+              !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
+              !assume that particles consist of ammonium sulfate at the surface (132g/mol) 
+              !and are dry at T = 20 deg. C and 50 perc. humidity
+              
+              !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
+              !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
+              sulfmmr_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &        ! [DRY part/kgA in bin it] 
+                   &  *(4./3.)*RPI*(mdw(it)/2.)**3.   &                   ! [mdw: dry diameter in m] 
+                   &  *dens_aer_dry                                       ! [dry aerosol mass density in kg/m3] 
+              
+              !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio) 
+              !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
+              sulfmmr(i,j) = sulfmmr(i,j) + sulfmmr_mode(i,j,it)
+              
+              !     nd_mode: particle concentration in different modes (DRY part/m3)
+              nd_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &             ! [DRY part/kgA in bin it] 
+                   & *pplay(i,j)/t_seri(i,j)/RD                           ! [air mass concentration in kg air /m3A]
+              
+              IF(flag_new_strat_compo) THEN
+                 !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
+                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)] 
+                      &  *4.*RPI*( mdw(it)*f_r_wetB(i,j,it)/2. )**2. &       ! [WET SA of part it in m2] 
+                      &  *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
+!    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
+                 samoment2 = samoment2 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
+                      &  *( mdw(it)*f_r_wetB(i,j,it)/2. )**2.                     ! [WET SA of part it in m2]
+!    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
+                 samoment3 = samoment3 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
+                      &  *( mdw(it)*f_r_wetB(i,j,it)/2. )**3.                     ! [WET SA of part it in m2]
+              ELSE
+!     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
+                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)] 
+                      &  *4.*RPI*( mdw(it)*f_r_wet(i,j)/2. )**2. &           ! [WET SA of part it in m2] 
+                      &  *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
+!    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
+                 samoment2 = samoment2 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
+                      &  *( mdw(it)*f_r_wet(i,j)/2. )**2.                          ! [WET SA of part it in m2]
+!    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
+                 samoment3 = samoment3 + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
+                      &  *( mdw(it)*f_r_wet(i,j)/2. )**3.                          ! [WET SA of part it in m2]
+              ENDIF
+           ENDDO
+!     reff_sulfate: effective radius of WET sulfate aerosols (cm)
+           reff_sulfate(i,j) = (samoment3 / samoment2) &
+                & *1.e2                                              ! conversion from m to cm
+        ENDDO
+     ENDDO
+     
   END SUBROUTINE traccoag
 
