Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 2155)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 2156)
@@ -30,5 +30,5 @@
       USE newsedim_mod, ONLY: newsedim
       USE datafile_mod, ONLY: datadir
-      USE improvedCO2clouds_mod, ONLY: improvedCO2clouds
+      USE improvedco2clouds_mod, ONLY: improvedco2clouds
 
       IMPLICIT NONE
@@ -610,5 +610,5 @@
 c      2.  Main call to the cloud schemes:
 c ==============================================================================
-        CALL improvedCO2clouds(ngrid,nlay,microtimestep,
+        CALL improvedco2clouds(ngrid,nlay,microtimestep,
      &     pplay,pplev,pteff,sum_subpdt,
      &     pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
Index: trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds_mod.F	(revision 2155)
+++ 	(revision )
@@ -1,670 +1,0 @@
-      MODULE improvedCO2clouds_mod
-
-      IMPLICIT NONE
-
-      CONTAINS
-
-      subroutine improvedCO2clouds(ngrid,nlay,microtimestep,
-     &             pplay,pplev,pteff,sum_subpdt,
-     &             pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
-     &             nq,tauscaling,
-     &             mem_Mccn_co2,mem_Mh2o_co2,mem_Nccn_co2,
-     &             No_dust,Mo_dust)
-      USE comcstfi_h, only: pi, g, cpp
-      USE updaterad, only: updaterice_micro, updaterice_microco2,
-     &                     updaterccnCO2
-      use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
-     &                      igcm_h2o_ice, igcm_ccn_mass,
-     &                      igcm_ccn_number, nuice_sed,
-     &                      igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
-     &                      igcm_ccnco2_number, nuiceco2_sed,
-     &                      rho_ice_co2,nuiceco2_ref
-      use conc_mod, only: mmean
-      use datafile_mod, only: datadir
-
-      implicit none
-      
-c------------------------------------------------------------------
-c  This routine is used to form CO2 clouds when a parcel of the GCM is
-c    saturated. It includes the ability to have supersaturation, a
-c    computation of the nucleation rates, growthrates and the
-c    scavenging of dust particles by clouds.
-c  It is worth noting that the amount of dust is computed using the
-c    dust optical depth computed in aeropacity.F. That's why
-c    the variable called "tauscaling" is used to convert
-c    pq(dust_mass) and pq(dust_number), which are relative
-c    quantities, to absolute and realistic quantities stored in zq.
-c    This has to be done to convert the inputs into absolute
-c    values, but also to convert the outputs back into relative
-c    values which are then used by the sedimentation and advection
-c    schemes.
-c CO2 ice particles can nucleate on both dust and on water ice particles
-c When CO2 ice is deposited onto a water ice particles, the particle is
-c removed from the water tracers.
-c Memory of the origin of the co2 particles is kept and thus the 
-c water cycle shouldn't be modified by this.
-c WARNING: no sedimentation of the water ice origin is performed
-c in the microphysical timestep in co2cloud.F. 
-
-c  Authors of the water ice clouds microphysics
-c J.-B. Madeleine, based on the work by Franck Montmessin
-c           (October 2011)
-c           T. Navarro, debug,correction, new scheme (October-April 2011)
-c           A. Spiga, optimization (February 2012)
-c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
-c of Constantino Listowski 
-c There is an energy limit to how much co2 can sublimate/condensate. It is 
-c defined by the difference of the GCM temperature with the co2 condensation 
-c temperature. 
-c Warning:
-c If meteoritic particles are activated and turn into co2 ice particles,
-c then they will be reversed in the dust tracers if the cloud sublimates
- 
-c------------------------------------------------------------------
-      include "callkeys.h"
-      include "microphys.h"
-c------------------------------------------------------------------
-c     Arguments:
-
-      INTEGER,INTENT(in) :: ngrid,nlay
-      integer,intent(in) :: nq         ! number of tracers
-      REAL,INTENT(in) :: microtimestep     ! physics time step (s)
-      REAL,INTENT(in) :: pplay(ngrid,nlay)     ! mid-layer pressure (Pa)
-      REAL,INTENT(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressure (Pa)
-      REAL,INTENT(in) :: pteff(ngrid,nlay) ! temperature at the middle of the
-                                 !   layers (K)
-      REAL,INTENT(in) :: sum_subpdt(ngrid,nlay) ! tendency on temperature from
-                                 !  previous physical parametrizations
-      REAL,INTENT(in) :: pqeff(ngrid,nlay,nq) ! tracers (kg/kg)
-      REAL,INTENT(in) :: sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers 
-                                 !  before condensation (kg/kg.s-1)
-      REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
-c     Outputs:
-      REAL,INTENT(out) :: subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
-                                   ! due to CO2 condensation (kg/kg.s-1)
-      ! condensation si igcm_co2_ice 
-      REAL,INTENT(out) :: subpdtcloudco2(ngrid,nlay)  ! tendency on temperature due
-                                   ! to latent heat
-      REAL,INTENT(out) :: No_dust(ngrid,nlay)
-      REAL,INTENT(out) :: Mo_dust(ngrid,nlay)
-
-c------------------------------------------------------------------
-c     Local variables:
-      LOGICAL,SAVE :: firstcall=.true.
-      REAL*8   derf ! Error function
-      INTEGER ig,l,i
-
-      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
-      REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
-                                ! used for nucleation of CO2 on ice-coated ccns
-      REAL rccnh2o(ngrid,nlay)    ! Water Ice mass mean radius (m)
-      REAL zq(ngrid,nlay,nq)  ! local value of tracers
-      REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
-      REAL zt(ngrid,nlay)       ! local value of temperature
-      REAL zqsat(ngrid,nlay)    ! saturation vapor pressure for CO2
-      real tcond(ngrid,nlay)
-      real zqco2(ngrid,nlay)
-      REAL lw                         !Latent heat of sublimation (J.kg-1) 
-      REAL,save :: l0,l1,l2,l3,l4
-      DOUBLE PRECISION dMice           ! mass of condensed ice
-      DOUBLE PRECISION sumcheck
-      DOUBLE PRECISION facteurmax!for energy limit on mass growth
-      DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
-      DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
-      DOUBLE PRECISION Mo ,No
-      
-c D.BARDET: sensibility test
-c      REAL, SAVE :: No ! when sensibility test
-c      DOUBLE PRECISION No_dust(ngrid,nlay) ! when sensibility test
-c      DOUBLE PRECISION Mo_dust(ngrid,nlay) ! when sensibility test 
-     
-      DOUBLE PRECISION  Rn, Rm, dev2,dev3, n_derf, m_derf
-      DOUBLE PRECISION mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
-      DOUBLE PRECISION mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
-      DOUBLE PRECISION mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
-      DOUBLE PRECISION interm1,interm2,interm3     
- 
-!     Radius used by the microphysical scheme (m)
-      DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
-      DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
-      DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
-
-      DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
-      DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
-      DOUBLE PRECISION rad_h2oice(nbinco2_cld) 
-
-c      REAL*8 sigco2      ! Co2-ice/air surface tension  (N.m)
-c      EXTERNAL sigco2
-
-      DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
-      DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
-
-      DOUBLE PRECISION rate(nbinco2_cld)  ! nucleation rate
-      DOUBLE PRECISION rateh2o(nbinco2_cld)  ! nucleation rate
-      REAL seq
-      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
-      DOUBLE PRECISION riceco2(ngrid,nlay)      ! CO2Ice mean radius (m)
-      REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
-                  
-      REAL rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
-      REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
-
-c      REAL res      ! Resistance growth
-      DOUBLE PRECISION Ic_rice      ! Mass transfer rate CO2 ice crystal
-      DOUBLE PRECISION ratioh2o_ccn
-      DOUBLE PRECISION vo2co2
-
-c     Parameters of the size discretization used by the microphysical scheme
-      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9   ! Minimum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6   ! Maximum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10  ! Minimum bounary radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4  ! Maximum boundary radius (m)
-      DOUBLE PRECISION vrat_cld                         ! Volume ratio
-      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1)         ! boundary values of each rad_cldco2 bin (m)
-      SAVE rb_cldco2
-      DOUBLE PRECISION dr_cld(nbinco2_cld)              ! width of each rad_cldco2 bin (m)
-      DOUBLE PRECISION vol_cld(nbinco2_cld)             ! particle volume for each bin (m3)
-
-      DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
-      REAL sigma_iceco2   ! Variance of the co2 ice and CCN distributions
-      SAVE sigma_iceco2 
-      REAL sigma_ice      ! Variance of the h2o ice and CCN distributions
-      SAVE sigma_ice
-      DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
-
-!     Variables for the meteoritic flux:
-      integer,parameter :: nbin_meteor=100
-      integer,parameter :: nlev_meteor=130
-      double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! 
-      double precision,save :: meteor(130,100)
-      double precision mtemp(100),pression_meteor(130)
-      logical file_ok
-      integer read_ok
-      integer nelem,lebon1,lebon2
-      double precision :: ltemp1(130),ltemp2(130)
-      integer ibin,j
-      integer,parameter :: uMeteor=666
-
-      IF (firstcall) THEN
-!=============================================================
-! 0. Definition of the size grid
-!=============================================================
-c       rad_cldco2 is the primary radius grid used for microphysics computation.
-c       The grid spacing is computed assuming a constant volume ratio
-c       between two consecutive bins; i.e. vrat_cld.
-c       vrat_cld is determined from the boundary values of the size grid: 
-c       rmin_cld and rmax_cld.
-c       The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
-c       dr_cld is the width of each rad_cldco2 bin.
-
-        vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
-        vrat_cld = exp(vrat_cld)
-        rb_cldco2(1)  = rbmin_cld
-        rad_cldco2(1) = rmin_cld
-        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
-        do i=1,nbinco2_cld-1
-          rad_cldco2(i+1)  = rad_cldco2(i) * vrat_cld**(1./3.)
-          vol_cld(i+1)  = vol_cld(i) * vrat_cld
-        enddo        
-        do i=1,nbinco2_cld
-          rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
-     &      rad_cldco2(i)
-          dr_cld(i)  = rb_cldco2(i+1) - rb_cldco2(i)
-        enddo
-        rb_cldco2(nbinco2_cld+1) = rbmax_cld
-        dr_cld(nbinco2_cld)   = rb_cldco2(nbinco2_cld+1) -
-     &       rb_cldco2(nbinco2_cld)
-        print*, ' '
-        print*,'Microphysics co2: size bin information:'
-        print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
-        print*,'-----------------------------------'
-        do i=1,nbinco2_cld
-          write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
-     &      dr_cld(i)
-        enddo
-        write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
-        print*,'-----------------------------------'
-        do i=1,nbinco2_cld+1
-            rb_cldco2(i) = log(rb_cldco2(i))  !! we save that so that it is not computed
-                                         !! at each timestep and gridpoint
-        enddo
-c       Contact parameter of co2 ice on dst ( m=cos(theta) )
-c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c        mteta  = 0.952
-c        mtetaco2 = 0.952
-c        write(*,*) 'co2_param contact parameter:', mtetaco2
-
-c       Volume of a co2 molecule (m3)
-        vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
-
-c       Variance of the ice and CCN distributions
-        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
-        sigma_ice = sqrt(log(1.+nuice_sed))
-
-  
-        write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
-        write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed
-        write(*,*) 'Volume of a co2 molecule:', vo1co2
-        
-        if (co2useh2o) then
-           write(*,*)
-           write(*,*) "co2useh2o=.true. in callphys.def"
-           write(*,*) "This means water ice particles can "
-           write(*,*) "serve as CCN for CO2 microphysics"
-        endif
-
-        if (meteo_flux) then
-           write(*,*)
-           write(*,*) "meteo_flux=.true. in callphys.def"
-           write(*,*) "meteoritic dust particles are available"
-           write(*,*) "for co2 ice nucleation! "
-           write(*,*) "Flux given by J. Plane (pressions,size bins)"
-           ! Initialisation of the flux: it is constant and is it saved
-           !We must interpolate the table to the GCM pressures
-           INQUIRE(FILE=TRIM(datadir)//
-     &       '/Meteo_flux_Plane.dat', EXIST=file_ok)
-           IF (.not. file_ok) THEN 
-              write(*,*) 'file Meteo_flux_Plane.dat should be in '
-     &             ,trim(datadir)
-              STOP
-           endif
-!used Variables
-           open(unit=uMeteor,file=trim(datadir)//
-     &          '/Meteo_flux_Plane.dat'
-     &          ,FORM='formatted')
-!13000 records (130 pressions x 100 bin sizes)
-           read(uMeteor,*) !skip 1 line
-           do i=1,130 
-              read(uMeteor,*,iostat=read_ok) pression_meteor(i)
-              if (read_ok==0) then
-                write(*,*) pression_meteor(i)
-              else
-                write(*,*) 'Error reading Meteo_flux_Plane.dat'
-                call abort_physic("CO2clouds",
-     &                    "Error reading Meteo_flux_Plane.dat",1)
-              endif
-           enddo
-           read(uMeteor,*)       !skip 1 line
-           do i=1,130 
-              do j=1,100        ! les mêmes 100 bins size que la distri nuclea: on touche pas
-                 read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
-                 if (read_ok/=0) then
-                   write(*,*) 'Error reading Meteo_flux_Plane.dat'
-                   call abort_physic("CO2clouds",
-     &                    "Error reading Meteo_flux_Plane.dat",1)
-                 endif
-              enddo
-!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
-           enddo
-           close(uMeteor)
-        write(*,*) "File meteo_flux read, end of firstcall in co2 micro"
-        endif                     !of if meteo_flux
-      
-        ! Parameter values for Latent heat computation
-        l0=595594d0      
-        l1=903.111d0    
-        l2=-11.5959d0   
-        l3=0.0528288d0
-        l4=-0.000103183d0
-        
-c D.BARDET:
-c        No = 1.e10        
-        firstcall=.false.
-      END IF
-!=============================================================
-! 1. Initialisation
-!=============================================================
-
-      meteor_ccn(:,:,:)=0.
-      rice(:,:) = 1.e-8
-      riceco2(:,:) = 1.e-11
-
-c     Initialize the tendencies
-      subpdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
-      subpdtcloudco2(1:ngrid,1:nlay)=0.
-      
-c pteff temperature layer; sum_subpdt dT.s-1
-c pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
-      zt(1:ngrid,1:nlay) = 
-     &      pteff(1:ngrid,1:nlay) + 
-     &      sum_subpdt(1:ngrid,1:nlay) * microtimestep 
-      zq(1:ngrid,1:nlay,1:nq) = 
-     &      pqeff(1:ngrid,1:nlay,1:nq) + 
-     &      sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
-      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
-     &     zq(1:ngrid,1:nlay,1:nq) = 1.e-30
-         
-         zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
-!=============================================================
-! 2. Compute saturation
-!=============================================================
-      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
-      dev3 = 1. / ( sqrt(2.) * sigma_ice )
-      call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
-      zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
-      CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
-!=============================================================
-! 3. Bonus: additional meteoritic particles for nucleation
-!=============================================================
-      if (meteo_flux) then
-         !pression_meteo(130)
-         !pplev(ngrid,nlay+1)
-         !meteo(130,100)
-         !resultat: meteo_ccn(ngrid,nlay,100)
-         do l=1,nlay
-            do ig=1,ngrid
-               masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
-               ltemp1=abs(pression_meteor(:)-pplev(ig,l))
-               ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
-               lebon1=minloc(ltemp1,DIM=1) 
-               lebon2=minloc(ltemp2,DIM=1)
-               nelem=lebon2-lebon1+1.
-               mtemp(:)=0d0     !mtemp(100) : valeurs pour les 100bins
-               do ibin=1,100
-                  mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
-               enddo
-               meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
-csi par m carre, x epaisseur/masse pour par kg/air 
-               !write(*,*) "masse air ig l=",masse(ig,l)
-               !check original unit with J. Plane
-            enddo
-         enddo
-      endif
-c ------------------------------------------------------------------------
-c ---------  Actual microphysics : Main loop over the GCM's grid ---------
-c ------------------------------------------------------------------------
-       DO l=1,nlay
-         DO ig=1,ngrid
-c       Get the partial pressure of co2 vapor and its saturation ratio
-           pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
-           satu = pco2 / zqsat(ig,l)
-
-           rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l)
-     &          -2.87e-6*zt(ig,l)*zt(ig,l)) !T-dependant CO2 ice density
-           vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
-           rho_ice_co2=rho_ice_co2T(ig,l)
-
-!=============================================================
-!4. Nucleation
-!=============================================================
-           IF ( satu .ge. 1 ) THEN ! if there is condensation
-
-          call updaterccnCO2(zq(ig,l,igcm_dust_mass),
-     &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
-
-c D.BARDET: sensibility test
-c              rdust=2.e-6
-      
-c Expand the dust moments into a binned distribution
-              
-              n_aer(:)=0.
-              m_aer(:)=0.
-
-              Mo =4.*pi*rho_dust*No*rdust(ig,l)**(3.)
-     &           *exp(9.*nuiceco2_ref/2.)/3. ! in Madeleine et al 2011
-     
-              No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
-
-              No_dust=No
-              Mo_dust=Mo
-
-              Rn = rdust(ig,l)
-              Rn = -log(Rn) 
-              Rm = Rn - 3. * sigma_iceco2*sigma_iceco2  
-              n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
-              m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
-
-              do i = 1, nbinco2_cld
-                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
-                 m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
-                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
-                 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
-                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
-                 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
-              enddo
-
-c Ajout meteor_ccn particles aux particules de poussière background
-              if (meteo_flux) then
-                 do i = 1, nbinco2_cld
-                    n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i) 
-                    m_aer(i) = m_aer(i) + 4./3.*pi*rho_dust
-     &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
-     &                *rad_cldco2(i)
-                 enddo
-              endif
-              
-c Same but with h2o particles as CCN only if co2useh2o=.true.
-              
-              n_aer_h2oice(:)=0.
-              m_aer_h2oice(:)=0.
-
-              if (co2useh2o) then
-                call updaterice_micro(zq(ig,l,igcm_h2o_ice),
-     &               zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
-     &                 tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-                Mo = zq(ig,l,igcm_h2o_ice) +
-     &               zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30 
-                     ! Total mass of H20 crystals,CCN included
-                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
-                Rn = rice(ig,l)
-                Rn = -log(Rn) 
-                Rm = Rn - 3. * sigma_ice*sigma_ice  
-                n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
-                m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
-                do i = 1, nbinco2_cld
-                  n_aer_h2oice(i) = -0.5 * No * n_derf 
-                  m_aer_h2oice(i) = -0.5 * Mo * m_derf 
-                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
-                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
-                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf 
-                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 
-                  rad_h2oice(i) = rad_cldco2(i)
-                enddo
-              endif
-             
-
-! Call to nucleation routine
-              call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
-     &             ,n_aer,rate,n_aer_h2oice
-     &             ,rad_h2oice,rateh2o,vo2co2)
-              dN = 0.
-              dM = 0.
-              dNh2o = 0.
-              dMh2o = 0.
-              do i = 1, nbinco2_cld
-                 Proba    =1.0-exp(-1.*microtimestep*rate(i))
-                 dN       = dN + n_aer(i) * Proba
-                 dM       = dM + m_aer(i) * Proba             
-              enddo
-              if (co2useh2o) then
-                 do i = 1, nbinco2_cld
-                    Probah2o = 1.0-exp(-1.*microtimestep*rateh2o(i))
-                    dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
-                    dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
-                 enddo
-              endif
-
-! dM  masse activée (kg) et dN nb particules par  kg d'air
-! Now increment CCN tracers and update dust tracers
-              dNN= dN/tauscaling(ig)
-              dMM= dM/tauscaling(ig)
-              dNN=min(dNN,zq(ig,l,igcm_dust_number))
-              dMM=min(dMM,zq(ig,l,igcm_dust_mass))
-              zq(ig,l,igcm_ccnco2_mass)   = 
-     &             zq(ig,l,igcm_ccnco2_mass)   + dMM
-              zq(ig,l,igcm_ccnco2_number) =
-     &             zq(ig,l,igcm_ccnco2_number) + dNN
-              zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM 
-              zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
-
-
-c Update CCN for CO2 nucleating on H2O CCN :
-              ! Warning: must keep memory of it 
-              if (co2useh2o) then
-                 dNNh2o=dNh2o/tauscaling(ig)
-                 dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
-                 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
-     &                +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))  
-                 dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn
-                 dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)*
-     &                tauscaling(ig)*ratioh2o_ccn
-                 dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
-                 dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
-                 dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
-                 zq(ig,l,igcm_ccnco2_mass)   = 
-     &                zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice+dMh2o_ccn
-                 zq(ig,l,igcm_ccnco2_number) = 
-     &                zq(ig,l,igcm_ccnco2_number) + dNNh2o
-                zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o 
-                zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
-                zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
-                mem_Mh2o_co2(ig,l)=mem_Mh2o_co2(ig,l)+dMh2o_ice
-                mem_Mccn_co2(ig,l)=mem_Mccn_co2(ig,l)+dMh2o_ccn
-                mem_Nccn_co2(ig,l)=mem_Nccn_co2(ig,l)+dNNh2o
-             endif ! of if co2useh2o
-           ENDIF   ! of is satu >1
-
-!=============================================================
-! 5. Ice growth: scheme for radius evolution
-!=============================================================
-
-c We trigger crystal growth if and only if there is at least one nuclei (N>1).
-c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
-c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
-             IF (zq(ig,l,igcm_ccnco2_number)
-     &         * tauscaling(ig)+1.e-30.ge. 1)THEN   ! we trigger crystal growth
-              
-              call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)),
-     &            dble(zq(ig,l,igcm_ccnco2_mass)),
-     &            dble(zq(ig,l,igcm_ccnco2_number)),
-     &            tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
-
-              Ic_rice=0.
-              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
-     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
-              facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
-              !specific heat of co2 ice = 1000 J.kg-1.K-1 
-              !specific heat of atm cpp = 744.5 J.kg-1.K-1
-
-c call scheme of microphys. mass growth for CO2
-              call massflowrateCO2(pplay(ig,l),zt(ig,l),
-     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice) 
-c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance ! 
-             
-              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then 
-                 Ic_rice=0.
-                 subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
-                 dMice=0
-                 
-              else
-                 dMice=zq(ig,l,igcm_ccnco2_number)*Ic_rice*microtimestep
-     &                *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
-                 !kg.s-1 par particule * nb particule par kg air*s
-                 ! = kg par kg air 
-              
-              dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
-              dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
-! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
-! latent heat release       >0 if growth i.e. if dMice >0
-              subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
-! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
-              !Now update tracers 
-              zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
-              zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
-              endif
-
-!=============================================================
-! 6. Dust cores releasing if no more co2 ice :
-!=============================================================
-
-              if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
-! On sublime tout
-                 if (co2useh2o) then
-                   if (mem_Mccn_co2(ig,l) .gt. 0) then
-                    zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
-     &                   +mem_Mccn_co2(ig,l)
-                   endif
-                   if (mem_Mh2o_co2(ig,l) .gt. 0) then
-                    zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
-     &                   +mem_Mh2o_co2(ig,l)
-                   endif
-                 
-                   if (mem_Nccn_co2(ig,l) .gt. 0) then
-                    zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
-     &                   +mem_Nccn_co2(ig,l)
-                   endif
-                 endif
-                    zq(ig,l,igcm_dust_mass) = 
-     &                   zq(ig,l,igcm_dust_mass)
-     &                   + zq(ig,l,igcm_ccnco2_mass)-
-     &                   (mem_Mh2o_co2(ig,l)+mem_Mccn_co2(ig,l))
-                    zq(ig,l,igcm_dust_number) = 
-     &                   zq(ig,l,igcm_dust_number)
-     &                   + zq(ig,l,igcm_ccnco2_number)
-     &                   -mem_Nccn_co2(ig,l)
-                 
-                    zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 
-     &                   + zq(ig,l,igcm_co2_ice)
-                 
-                 zq(ig,l,igcm_ccnco2_mass)=0.
-                 zq(ig,l,igcm_co2_ice)=0.
-                 zq(ig,l,igcm_ccnco2_number)=0.
-                 mem_Nccn_co2(ig,l)=0.
-                 mem_Mh2o_co2(ig,l)=0.
-                 mem_Mccn_co2(ig,l)=0.
-                 riceco2(ig,l)=0.
-
-              endif !of if co2_ice <1e-25
-
-              ENDIF            ! of if NCCN > 1
-          ENDDO                ! of ig loop
-        ENDDO                  ! of nlayer loop  
-
-!=============================================================
-! 7. END: get cloud tendencies 
-!=============================================================
-
-          ! Get cloud tendencies
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
-     &       (zq(1:ngrid,1:nlay,igcm_co2) - 
-     &       zq0(1:ngrid,1:nlay,igcm_co2))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
-     &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
-     &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
-        
-        if (co2useh2o) then
-           subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
-     &         (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
-     &         zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
-           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
-     &         (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
-     &         zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
-           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
-     &         (zq(1:ngrid,1:nlay,igcm_ccn_number) -
-     &         zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
-        endif
-
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
-     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
-     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
-     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
-     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/microtimestep
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
-     &       (zq(1:ngrid,1:nlay,igcm_dust_mass) -
-     &       zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep 
-        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
-     &       (zq(1:ngrid,1:nlay,igcm_dust_number) -
-     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
-
-
-c     TEST D.BARDET 
-      call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussière"
-     &        ,"part/kg",3,No_dust)
-      call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussière"
-     &        ,"kg/kg ",3,Mo_dust)      
-
-        END SUBROUTINE improvedCO2clouds
-
-        END MODULE improvedCO2clouds_mod
-
Index: trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F	(revision 2156)
+++ trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F	(revision 2156)
@@ -0,0 +1,670 @@
+      MODULE improvedco2clouds_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      subroutine improvedco2clouds(ngrid,nlay,microtimestep,
+     &             pplay,pplev,pteff,sum_subpdt,
+     &             pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
+     &             nq,tauscaling,
+     &             mem_Mccn_co2,mem_Mh2o_co2,mem_Nccn_co2,
+     &             No_dust,Mo_dust)
+      USE comcstfi_h, only: pi, g, cpp
+      USE updaterad, only: updaterice_micro, updaterice_microco2,
+     &                     updaterccnCO2
+      use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
+     &                      igcm_h2o_ice, igcm_ccn_mass,
+     &                      igcm_ccn_number, nuice_sed,
+     &                      igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
+     &                      igcm_ccnco2_number, nuiceco2_sed,
+     &                      rho_ice_co2,nuiceco2_ref
+      use conc_mod, only: mmean
+      use datafile_mod, only: datadir
+
+      implicit none
+      
+c------------------------------------------------------------------
+c  This routine is used to form CO2 clouds when a parcel of the GCM is
+c    saturated. It includes the ability to have supersaturation, a
+c    computation of the nucleation rates, growthrates and the
+c    scavenging of dust particles by clouds.
+c  It is worth noting that the amount of dust is computed using the
+c    dust optical depth computed in aeropacity.F. That's why
+c    the variable called "tauscaling" is used to convert
+c    pq(dust_mass) and pq(dust_number), which are relative
+c    quantities, to absolute and realistic quantities stored in zq.
+c    This has to be done to convert the inputs into absolute
+c    values, but also to convert the outputs back into relative
+c    values which are then used by the sedimentation and advection
+c    schemes.
+c CO2 ice particles can nucleate on both dust and on water ice particles
+c When CO2 ice is deposited onto a water ice particles, the particle is
+c removed from the water tracers.
+c Memory of the origin of the co2 particles is kept and thus the 
+c water cycle shouldn't be modified by this.
+c WARNING: no sedimentation of the water ice origin is performed
+c in the microphysical timestep in co2cloud.F. 
+
+c  Authors of the water ice clouds microphysics
+c J.-B. Madeleine, based on the work by Franck Montmessin
+c           (October 2011)
+c           T. Navarro, debug,correction, new scheme (October-April 2011)
+c           A. Spiga, optimization (February 2012)
+c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
+c of Constantino Listowski 
+c There is an energy limit to how much co2 can sublimate/condensate. It is 
+c defined by the difference of the GCM temperature with the co2 condensation 
+c temperature. 
+c Warning:
+c If meteoritic particles are activated and turn into co2 ice particles,
+c then they will be reversed in the dust tracers if the cloud sublimates
+ 
+c------------------------------------------------------------------
+      include "callkeys.h"
+      include "microphys.h"
+c------------------------------------------------------------------
+c     Arguments:
+
+      INTEGER,INTENT(in) :: ngrid,nlay
+      integer,intent(in) :: nq         ! number of tracers
+      REAL,INTENT(in) :: microtimestep     ! physics time step (s)
+      REAL,INTENT(in) :: pplay(ngrid,nlay)     ! mid-layer pressure (Pa)
+      REAL,INTENT(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressure (Pa)
+      REAL,INTENT(in) :: pteff(ngrid,nlay) ! temperature at the middle of the
+                                 !   layers (K)
+      REAL,INTENT(in) :: sum_subpdt(ngrid,nlay) ! tendency on temperature from
+                                 !  previous physical parametrizations
+      REAL,INTENT(in) :: pqeff(ngrid,nlay,nq) ! tracers (kg/kg)
+      REAL,INTENT(in) :: sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers 
+                                 !  before condensation (kg/kg.s-1)
+      REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
+c     Outputs:
+      REAL,INTENT(out) :: subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
+                                   ! due to CO2 condensation (kg/kg.s-1)
+      ! condensation si igcm_co2_ice 
+      REAL,INTENT(out) :: subpdtcloudco2(ngrid,nlay)  ! tendency on temperature due
+                                   ! to latent heat
+      REAL,INTENT(out) :: No_dust(ngrid,nlay)
+      REAL,INTENT(out) :: Mo_dust(ngrid,nlay)
+
+c------------------------------------------------------------------
+c     Local variables:
+      LOGICAL,SAVE :: firstcall=.true.
+      REAL*8   derf ! Error function
+      INTEGER ig,l,i
+
+      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
+      REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
+                                ! used for nucleation of CO2 on ice-coated ccns
+      REAL rccnh2o(ngrid,nlay)    ! Water Ice mass mean radius (m)
+      REAL zq(ngrid,nlay,nq)  ! local value of tracers
+      REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
+      REAL zt(ngrid,nlay)       ! local value of temperature
+      REAL zqsat(ngrid,nlay)    ! saturation vapor pressure for CO2
+      real tcond(ngrid,nlay)
+      real zqco2(ngrid,nlay)
+      REAL lw                         !Latent heat of sublimation (J.kg-1) 
+      REAL,save :: l0,l1,l2,l3,l4
+      DOUBLE PRECISION dMice           ! mass of condensed ice
+      DOUBLE PRECISION sumcheck
+      DOUBLE PRECISION facteurmax!for energy limit on mass growth
+      DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
+      DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
+      DOUBLE PRECISION Mo ,No
+      
+c D.BARDET: sensibility test
+c      REAL, SAVE :: No ! when sensibility test
+c      DOUBLE PRECISION No_dust(ngrid,nlay) ! when sensibility test
+c      DOUBLE PRECISION Mo_dust(ngrid,nlay) ! when sensibility test 
+     
+      DOUBLE PRECISION  Rn, Rm, dev2,dev3, n_derf, m_derf
+      DOUBLE PRECISION mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
+      DOUBLE PRECISION mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
+      DOUBLE PRECISION mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
+      DOUBLE PRECISION interm1,interm2,interm3     
+ 
+!     Radius used by the microphysical scheme (m)
+      DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
+      DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
+      DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
+
+      DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
+      DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
+      DOUBLE PRECISION rad_h2oice(nbinco2_cld) 
+
+c      REAL*8 sigco2      ! Co2-ice/air surface tension  (N.m)
+c      EXTERNAL sigco2
+
+      DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
+      DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
+
+      DOUBLE PRECISION rate(nbinco2_cld)  ! nucleation rate
+      DOUBLE PRECISION rateh2o(nbinco2_cld)  ! nucleation rate
+      REAL seq
+      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
+      DOUBLE PRECISION riceco2(ngrid,nlay)      ! CO2Ice mean radius (m)
+      REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
+                  
+      REAL rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
+      REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
+
+c      REAL res      ! Resistance growth
+      DOUBLE PRECISION Ic_rice      ! Mass transfer rate CO2 ice crystal
+      DOUBLE PRECISION ratioh2o_ccn
+      DOUBLE PRECISION vo2co2
+
+c     Parameters of the size discretization used by the microphysical scheme
+      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9   ! Minimum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6   ! Maximum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10  ! Minimum bounary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4  ! Maximum boundary radius (m)
+      DOUBLE PRECISION vrat_cld                         ! Volume ratio
+      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1)         ! boundary values of each rad_cldco2 bin (m)
+      SAVE rb_cldco2
+      DOUBLE PRECISION dr_cld(nbinco2_cld)              ! width of each rad_cldco2 bin (m)
+      DOUBLE PRECISION vol_cld(nbinco2_cld)             ! particle volume for each bin (m3)
+
+      DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
+      REAL sigma_iceco2   ! Variance of the co2 ice and CCN distributions
+      SAVE sigma_iceco2 
+      REAL sigma_ice      ! Variance of the h2o ice and CCN distributions
+      SAVE sigma_ice
+      DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
+
+!     Variables for the meteoritic flux:
+      integer,parameter :: nbin_meteor=100
+      integer,parameter :: nlev_meteor=130
+      double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! 
+      double precision,save :: meteor(130,100)
+      double precision mtemp(100),pression_meteor(130)
+      logical file_ok
+      integer read_ok
+      integer nelem,lebon1,lebon2
+      double precision :: ltemp1(130),ltemp2(130)
+      integer ibin,j
+      integer,parameter :: uMeteor=666
+
+      IF (firstcall) THEN
+!=============================================================
+! 0. Definition of the size grid
+!=============================================================
+c       rad_cldco2 is the primary radius grid used for microphysics computation.
+c       The grid spacing is computed assuming a constant volume ratio
+c       between two consecutive bins; i.e. vrat_cld.
+c       vrat_cld is determined from the boundary values of the size grid: 
+c       rmin_cld and rmax_cld.
+c       The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
+c       dr_cld is the width of each rad_cldco2 bin.
+
+        vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
+        vrat_cld = exp(vrat_cld)
+        rb_cldco2(1)  = rbmin_cld
+        rad_cldco2(1) = rmin_cld
+        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
+        do i=1,nbinco2_cld-1
+          rad_cldco2(i+1)  = rad_cldco2(i) * vrat_cld**(1./3.)
+          vol_cld(i+1)  = vol_cld(i) * vrat_cld
+        enddo        
+        do i=1,nbinco2_cld
+          rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
+     &      rad_cldco2(i)
+          dr_cld(i)  = rb_cldco2(i+1) - rb_cldco2(i)
+        enddo
+        rb_cldco2(nbinco2_cld+1) = rbmax_cld
+        dr_cld(nbinco2_cld)   = rb_cldco2(nbinco2_cld+1) -
+     &       rb_cldco2(nbinco2_cld)
+        print*, ' '
+        print*,'Microphysics co2: size bin information:'
+        print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
+        print*,'-----------------------------------'
+        do i=1,nbinco2_cld
+          write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
+     &      dr_cld(i)
+        enddo
+        write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
+        print*,'-----------------------------------'
+        do i=1,nbinco2_cld+1
+            rb_cldco2(i) = log(rb_cldco2(i))  !! we save that so that it is not computed
+                                         !! at each timestep and gridpoint
+        enddo
+c       Contact parameter of co2 ice on dst ( m=cos(theta) )
+c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c        mteta  = 0.952
+c        mtetaco2 = 0.952
+c        write(*,*) 'co2_param contact parameter:', mtetaco2
+
+c       Volume of a co2 molecule (m3)
+        vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
+
+c       Variance of the ice and CCN distributions
+        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
+        sigma_ice = sqrt(log(1.+nuice_sed))
+
+  
+        write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
+        write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed
+        write(*,*) 'Volume of a co2 molecule:', vo1co2
+        
+        if (co2useh2o) then
+           write(*,*)
+           write(*,*) "co2useh2o=.true. in callphys.def"
+           write(*,*) "This means water ice particles can "
+           write(*,*) "serve as CCN for CO2 microphysics"
+        endif
+
+        if (meteo_flux) then
+           write(*,*)
+           write(*,*) "meteo_flux=.true. in callphys.def"
+           write(*,*) "meteoritic dust particles are available"
+           write(*,*) "for co2 ice nucleation! "
+           write(*,*) "Flux given by J. Plane (pressions,size bins)"
+           ! Initialisation of the flux: it is constant and is it saved
+           !We must interpolate the table to the GCM pressures
+           INQUIRE(FILE=TRIM(datadir)//
+     &       '/Meteo_flux_Plane.dat', EXIST=file_ok)
+           IF (.not. file_ok) THEN 
+              write(*,*) 'file Meteo_flux_Plane.dat should be in '
+     &             ,trim(datadir)
+              STOP
+           endif
+!used Variables
+           open(unit=uMeteor,file=trim(datadir)//
+     &          '/Meteo_flux_Plane.dat'
+     &          ,FORM='formatted')
+!13000 records (130 pressions x 100 bin sizes)
+           read(uMeteor,*) !skip 1 line
+           do i=1,130 
+              read(uMeteor,*,iostat=read_ok) pression_meteor(i)
+              if (read_ok==0) then
+                write(*,*) pression_meteor(i)
+              else
+                write(*,*) 'Error reading Meteo_flux_Plane.dat'
+                call abort_physic("CO2clouds",
+     &                    "Error reading Meteo_flux_Plane.dat",1)
+              endif
+           enddo
+           read(uMeteor,*)       !skip 1 line
+           do i=1,130 
+              do j=1,100        ! les mêmes 100 bins size que la distri nuclea: on touche pas
+                 read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
+                 if (read_ok/=0) then
+                   write(*,*) 'Error reading Meteo_flux_Plane.dat'
+                   call abort_physic("CO2clouds",
+     &                    "Error reading Meteo_flux_Plane.dat",1)
+                 endif
+              enddo
+!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
+           enddo
+           close(uMeteor)
+        write(*,*) "File meteo_flux read, end of firstcall in co2 micro"
+        endif                     !of if meteo_flux
+      
+        ! Parameter values for Latent heat computation
+        l0=595594d0      
+        l1=903.111d0    
+        l2=-11.5959d0   
+        l3=0.0528288d0
+        l4=-0.000103183d0
+        
+c D.BARDET:
+c        No = 1.e10        
+        firstcall=.false.
+      END IF
+!=============================================================
+! 1. Initialisation
+!=============================================================
+
+      meteor_ccn(:,:,:)=0.
+      rice(:,:) = 1.e-8
+      riceco2(:,:) = 1.e-11
+
+c     Initialize the tendencies
+      subpdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
+      subpdtcloudco2(1:ngrid,1:nlay)=0.
+      
+c pteff temperature layer; sum_subpdt dT.s-1
+c pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
+      zt(1:ngrid,1:nlay) = 
+     &      pteff(1:ngrid,1:nlay) + 
+     &      sum_subpdt(1:ngrid,1:nlay) * microtimestep 
+      zq(1:ngrid,1:nlay,1:nq) = 
+     &      pqeff(1:ngrid,1:nlay,1:nq) + 
+     &      sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
+      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
+     &     zq(1:ngrid,1:nlay,1:nq) = 1.e-30
+         
+         zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
+!=============================================================
+! 2. Compute saturation
+!=============================================================
+      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
+      dev3 = 1. / ( sqrt(2.) * sigma_ice )
+      call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
+      zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
+      CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
+!=============================================================
+! 3. Bonus: additional meteoritic particles for nucleation
+!=============================================================
+      if (meteo_flux) then
+         !pression_meteo(130)
+         !pplev(ngrid,nlay+1)
+         !meteo(130,100)
+         !resultat: meteo_ccn(ngrid,nlay,100)
+         do l=1,nlay
+            do ig=1,ngrid
+               masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+               ltemp1=abs(pression_meteor(:)-pplev(ig,l))
+               ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
+               lebon1=minloc(ltemp1,DIM=1) 
+               lebon2=minloc(ltemp2,DIM=1)
+               nelem=lebon2-lebon1+1.
+               mtemp(:)=0d0     !mtemp(100) : valeurs pour les 100bins
+               do ibin=1,100
+                  mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
+               enddo
+               meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
+csi par m carre, x epaisseur/masse pour par kg/air 
+               !write(*,*) "masse air ig l=",masse(ig,l)
+               !check original unit with J. Plane
+            enddo
+         enddo
+      endif
+c ------------------------------------------------------------------------
+c ---------  Actual microphysics : Main loop over the GCM's grid ---------
+c ------------------------------------------------------------------------
+       DO l=1,nlay
+         DO ig=1,ngrid
+c       Get the partial pressure of co2 vapor and its saturation ratio
+           pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
+           satu = pco2 / zqsat(ig,l)
+
+           rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l)
+     &          -2.87e-6*zt(ig,l)*zt(ig,l)) !T-dependant CO2 ice density
+           vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
+           rho_ice_co2=rho_ice_co2T(ig,l)
+
+!=============================================================
+!4. Nucleation
+!=============================================================
+           IF ( satu .ge. 1 ) THEN ! if there is condensation
+
+          call updaterccnCO2(zq(ig,l,igcm_dust_mass),
+     &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
+
+c D.BARDET: sensibility test
+c              rdust=2.e-6
+      
+c Expand the dust moments into a binned distribution
+              
+              n_aer(:)=0.
+              m_aer(:)=0.
+
+              Mo =4.*pi*rho_dust*No*rdust(ig,l)**(3.)
+     &           *exp(9.*nuiceco2_ref/2.)/3. ! in Madeleine et al 2011
+     
+              No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
+
+              No_dust=No
+              Mo_dust=Mo
+
+              Rn = rdust(ig,l)
+              Rn = -log(Rn) 
+              Rm = Rn - 3. * sigma_iceco2*sigma_iceco2  
+              n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
+              m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
+
+              do i = 1, nbinco2_cld
+                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
+                 m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
+                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
+                 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
+                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
+                 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
+              enddo
+
+c Ajout meteor_ccn particles aux particules de poussière background
+              if (meteo_flux) then
+                 do i = 1, nbinco2_cld
+                    n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i) 
+                    m_aer(i) = m_aer(i) + 4./3.*pi*rho_dust
+     &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
+     &                *rad_cldco2(i)
+                 enddo
+              endif
+              
+c Same but with h2o particles as CCN only if co2useh2o=.true.
+              
+              n_aer_h2oice(:)=0.
+              m_aer_h2oice(:)=0.
+
+              if (co2useh2o) then
+                call updaterice_micro(zq(ig,l,igcm_h2o_ice),
+     &               zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
+     &                 tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+                Mo = zq(ig,l,igcm_h2o_ice) +
+     &               zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30 
+                     ! Total mass of H20 crystals,CCN included
+                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
+                Rn = rice(ig,l)
+                Rn = -log(Rn) 
+                Rm = Rn - 3. * sigma_ice*sigma_ice  
+                n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
+                m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
+                do i = 1, nbinco2_cld
+                  n_aer_h2oice(i) = -0.5 * No * n_derf 
+                  m_aer_h2oice(i) = -0.5 * Mo * m_derf 
+                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
+                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
+                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf 
+                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 
+                  rad_h2oice(i) = rad_cldco2(i)
+                enddo
+              endif
+             
+
+! Call to nucleation routine
+              call nucleaco2(dble(pco2),zt(ig,l),dble(satu)
+     &             ,n_aer,rate,n_aer_h2oice
+     &             ,rad_h2oice,rateh2o,vo2co2)
+              dN = 0.
+              dM = 0.
+              dNh2o = 0.
+              dMh2o = 0.
+              do i = 1, nbinco2_cld
+                 Proba    =1.0-exp(-1.*microtimestep*rate(i))
+                 dN       = dN + n_aer(i) * Proba
+                 dM       = dM + m_aer(i) * Proba             
+              enddo
+              if (co2useh2o) then
+                 do i = 1, nbinco2_cld
+                    Probah2o = 1.0-exp(-1.*microtimestep*rateh2o(i))
+                    dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
+                    dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
+                 enddo
+              endif
+
+! dM  masse activée (kg) et dN nb particules par  kg d'air
+! Now increment CCN tracers and update dust tracers
+              dNN= dN/tauscaling(ig)
+              dMM= dM/tauscaling(ig)
+              dNN=min(dNN,zq(ig,l,igcm_dust_number))
+              dMM=min(dMM,zq(ig,l,igcm_dust_mass))
+              zq(ig,l,igcm_ccnco2_mass)   = 
+     &             zq(ig,l,igcm_ccnco2_mass)   + dMM
+              zq(ig,l,igcm_ccnco2_number) =
+     &             zq(ig,l,igcm_ccnco2_number) + dNN
+              zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM 
+              zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
+
+
+c Update CCN for CO2 nucleating on H2O CCN :
+              ! Warning: must keep memory of it 
+              if (co2useh2o) then
+                 dNNh2o=dNh2o/tauscaling(ig)
+                 dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
+                 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
+     &                +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))  
+                 dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn
+                 dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)*
+     &                tauscaling(ig)*ratioh2o_ccn
+                 dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
+                 dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
+                 dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
+                 zq(ig,l,igcm_ccnco2_mass)   = 
+     &                zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice+dMh2o_ccn
+                 zq(ig,l,igcm_ccnco2_number) = 
+     &                zq(ig,l,igcm_ccnco2_number) + dNNh2o
+                zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o 
+                zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
+                zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
+                mem_Mh2o_co2(ig,l)=mem_Mh2o_co2(ig,l)+dMh2o_ice
+                mem_Mccn_co2(ig,l)=mem_Mccn_co2(ig,l)+dMh2o_ccn
+                mem_Nccn_co2(ig,l)=mem_Nccn_co2(ig,l)+dNNh2o
+             endif ! of if co2useh2o
+           ENDIF   ! of is satu >1
+
+!=============================================================
+! 5. Ice growth: scheme for radius evolution
+!=============================================================
+
+c We trigger crystal growth if and only if there is at least one nuclei (N>1).
+c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
+c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
+             IF (zq(ig,l,igcm_ccnco2_number)
+     &         * tauscaling(ig)+1.e-30.ge. 1)THEN   ! we trigger crystal growth
+              
+              call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)),
+     &            dble(zq(ig,l,igcm_ccnco2_mass)),
+     &            dble(zq(ig,l,igcm_ccnco2_number)),
+     &            tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
+
+              Ic_rice=0.
+              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
+     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
+              facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
+              !specific heat of co2 ice = 1000 J.kg-1.K-1 
+              !specific heat of atm cpp = 744.5 J.kg-1.K-1
+
+c call scheme of microphys. mass growth for CO2
+              call massflowrateco2(pplay(ig,l),zt(ig,l),
+     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice) 
+c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance ! 
+             
+              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then 
+                 Ic_rice=0.
+                 subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
+                 dMice=0
+                 
+              else
+                 dMice=zq(ig,l,igcm_ccnco2_number)*Ic_rice*microtimestep
+     &                *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
+                 !kg.s-1 par particule * nb particule par kg air*s
+                 ! = kg par kg air 
+              
+              dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
+              dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
+! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
+! latent heat release       >0 if growth i.e. if dMice >0
+              subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
+! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
+              !Now update tracers 
+              zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
+              zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
+              endif
+
+!=============================================================
+! 6. Dust cores releasing if no more co2 ice :
+!=============================================================
+
+              if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
+! On sublime tout
+                 if (co2useh2o) then
+                   if (mem_Mccn_co2(ig,l) .gt. 0) then
+                    zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
+     &                   +mem_Mccn_co2(ig,l)
+                   endif
+                   if (mem_Mh2o_co2(ig,l) .gt. 0) then
+                    zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
+     &                   +mem_Mh2o_co2(ig,l)
+                   endif
+                 
+                   if (mem_Nccn_co2(ig,l) .gt. 0) then
+                    zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
+     &                   +mem_Nccn_co2(ig,l)
+                   endif
+                 endif
+                    zq(ig,l,igcm_dust_mass) = 
+     &                   zq(ig,l,igcm_dust_mass)
+     &                   + zq(ig,l,igcm_ccnco2_mass)-
+     &                   (mem_Mh2o_co2(ig,l)+mem_Mccn_co2(ig,l))
+                    zq(ig,l,igcm_dust_number) = 
+     &                   zq(ig,l,igcm_dust_number)
+     &                   + zq(ig,l,igcm_ccnco2_number)
+     &                   -mem_Nccn_co2(ig,l)
+                 
+                    zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 
+     &                   + zq(ig,l,igcm_co2_ice)
+                 
+                 zq(ig,l,igcm_ccnco2_mass)=0.
+                 zq(ig,l,igcm_co2_ice)=0.
+                 zq(ig,l,igcm_ccnco2_number)=0.
+                 mem_Nccn_co2(ig,l)=0.
+                 mem_Mh2o_co2(ig,l)=0.
+                 mem_Mccn_co2(ig,l)=0.
+                 riceco2(ig,l)=0.
+
+              endif !of if co2_ice <1e-25
+
+              ENDIF            ! of if NCCN > 1
+          ENDDO                ! of ig loop
+        ENDDO                  ! of nlayer loop  
+
+!=============================================================
+! 7. END: get cloud tendencies 
+!=============================================================
+
+          ! Get cloud tendencies
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
+     &       (zq(1:ngrid,1:nlay,igcm_co2) - 
+     &       zq0(1:ngrid,1:nlay,igcm_co2))/microtimestep
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
+     &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
+     &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
+        
+        if (co2useh2o) then
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
+     &         (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
+     &         zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
+     &         (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
+     &         zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
+           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
+     &         (zq(1:ngrid,1:nlay,igcm_ccn_number) -
+     &         zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
+        endif
+
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/microtimestep
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/microtimestep
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
+     &       (zq(1:ngrid,1:nlay,igcm_dust_mass) -
+     &       zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep 
+        subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
+     &       (zq(1:ngrid,1:nlay,igcm_dust_number) -
+     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
+
+
+c     TEST D.BARDET 
+      call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussière"
+     &        ,"part/kg",3,No_dust)
+      call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussière"
+     &        ,"kg/kg ",3,Mo_dust)      
+
+        END SUBROUTINE improvedco2clouds
+
+        END MODULE improvedco2clouds_mod
+
Index: trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F	(revision 2155)
+++ 	(revision )
@@ -1,540 +1,0 @@
-c=======================================================================
-      subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
-c
-c     Determination of the mass transfer rate for CO2 condensation & 
-c     sublimation
-c
-c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
-c           particle radius (Radius), molecular mass of the atm (Matm)
-c   output:  MASS FLUX Ic
-c
-c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
-c  
-c
-c   Updates:
-c   --------
-c   December 2017 - C. Listowski - Simplification of the derivation of
-c   massflowrate by using explicit formula for surface temperature,
-c   No Newton-Raphson routine anymore- see comment at relevant line
-c=======================================================================
-      USE comcstfi_h, ONLY: pi
-
-      implicit none
-
-      include "microphys.h"
-
-c   arguments: INPUT
-c   ----------
-      REAL,INTENT(in) :: T,Matm
-      REAL*8,INTENT(in) :: SAT
-      REAL,INTENT(in) :: P
-      DOUBLE PRECISION,INTENT(in) :: Radius
-c   arguments: OUTPUT
-c   ----------
-      DOUBLE PRECISION,INTENT(out) ::   Ic
-c   Local Variables
-c   ----------
-      DOUBLE PRECISION   Tsurf
-      DOUBLE PRECISION   C0,C1,C2
-      DOUBLE PRECISION   kmix,Lsub,cond
-      DOUBLE PRECISION   Ak
-
-      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak)
-   
-      Tsurf = 1./C1*dlog(Sat/Ak) + T 
-
-      !Note by CL - dec 2017 (see also technical note)
-      !The above is a simplified version of Tsurf
-      !compared to the one used by Listowski et al. 2014 (Ta), where a
-      !Newton-Raphson routine must be used. Approximations made by
-      !considering the orders of magnitude of the different factors lead to
-      !simplification of the equation 5 of Listowski et al. (2014).
-      !The error compared to the exact value determined by NR iterations
-      !is less than 0.6% for all sizes, pressures, supersaturations
-      !relevant to present Mars. Should also be ok for most conditions
-      !in ancient Mars (However, needs to be double-cheked, as explained in
-      !(Listowski et al. 2013,JGR)
-
-      cond = 4.*pi*Radius*kmix
-      Ic = cond*(Tsurf-T)/Lsub
-  
-      END
-
-c********************************************************************************
-
-      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak)
-
-c********************************************************************************
-c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
-      use tracer_mod, only: rho_ice_co2
-      USE comcstfi_h, ONLY: pi
-
-      implicit none
-      include "microphys.h"
-
-c   arguments: INPUT
-c   ----------------
-      REAL,INTENT(in) :: P
-      REAL,INTENT(in) :: T
-      REAL*8,INTENT(in) :: S
-      DOUBLE PRECISION,INTENT(in) :: rc
-      REAL,INTENT(in) :: Matm !g.mol-1 ( = mmean(ig,l) )
-c   arguments: OUTPUT
-c   ----------
-      DOUBLE PRECISION,INTENT(out) ::  C0,C1,C2
-      DOUBLE PRECISION,INTENT(out) ::  kmix,Lsub
-
-c   local:
-c   ------
-      DOUBLE PRECISION Cpatm,Cpn2,Cpco2
-      DOUBLE PRECISION psat, xinf, pco2
-      DOUBLE PRECISION Dv     
-      DOUBLE PRECISION l0,l1,l2,l3,l4            
-      DOUBLE PRECISION knudsen, a, lambda      ! F and S correction
-      DOUBLE PRECISION Ak                       ! kelvin factor    
-      DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th
-
-
-c     DEFINE heat cap. J.kg-1.K-1 and To
-
-      data Cpco2/0.7e3/
-      data Cpn2/1e3/
-
-      kmix = 0d0
-      Lsub = 0d0
-
-      C0 = 0d0
-      C1 = 0d0
-      C2 = 0d0
-
-c     Equilibirum pressure over a flat surface
-      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
-c     Compute transport coefficient
-      pco2 = psat * dble(S)
-c     Latent heat of sublimation if CO2  co2 (J.kg-1)
-c     version Azreg_Ainou (J/kg) :
-      l0=595594.      
-      l1=903.111     
-      l2=-11.5959    
-      l3=0.0528288 
-      l4=-0.000103183
-      Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * 
-     &     dble(T)**3 + l4 * dble(T)**4 ! J/kg
-c     atmospheric density
-      rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
-      rhoatm = rhoatm * 1.00e-3 !kg.m-3
-      call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
-      call  Diffcoeff(P, T, Dv)  
-
-      Dv = Dv * 1.00e-4         !!! cm2.s-1  to m2.s-1
-
-c     ----- FS correction for Diff
-      vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1
-      knudsen = 3*Dv / (vthco2 * rc) 
-      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
-      Dv      = Dv / (1. + lambda * knudsen)
-c     ----- FS correction for Kth 
-      vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg 
-      Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1
-      lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
-     &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer
-      knudsen = lpmt / rc
-      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
-      kmix    = kmix /  (1. + lambda * knudsen)
-c     --------------------- ASSIGN coeff values for FUNCTION
-      xinf = dble(S) * psat / dble(P)
-      Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
-      C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
-     &     (rgp*dble(T)))
-      C1 = Lsub*mco2/(rgp*dble(T)**2)
-      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
-
-      END
-
-
-c======================================================================
-      subroutine Diffcoeff(P, T, Diff)
-c     Compute diffusion coefficient CO2/N2
-c     cited in Ilona's lecture - from Reid et al. 1987
-c======================================================================
-       IMPLICIT NONE
-
-       include "microphys.h"
-
-c      arguments
-c     -----------
-      
-      REAL,INTENT(in) :: P
-      REAL,INTENT(in) :: T
-      
-c     output
-c     -----------
-
-      DOUBLE PRECISION,INTENT(out) :: Diff
-    
-c      local
-c     -----------
-
-      REAL Pbar                     !!! has to be in bar for the formula
-      DOUBLE PRECISION dva, dvb, Mab  ! Mab has to be in g.mol-1
-        
-        Pbar = P * 1d-5
-    
-        Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000.
-
-  	dva = 26.9        ! diffusion volume of CO2,  Reid et al. 1987 (cited in Ilona's lecture)
-  	dvb = 18.5        ! diffusion volume of N2
-    
-        Diff  = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) 
-     &       * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1
-  
-       RETURN
-
-       END
-
-
-c======================================================================
-
-         subroutine KthMixNEW(Kthmix,T,x,rho)
-
-c        Compute thermal conductivity of CO2/N2 mixture
-c         (***WITHOUT*** USE OF VISCOSITY)
-
-c          (Mason & Saxena, 1958 - Wassiljeva 1904)
-c======================================================================
-
-       implicit none
-       
-       include "microphys.h"
-c      arguments
-c     -----------
-         
-         REAL,INTENT(in) :: T
-         DOUBLE PRECISION,INTENT(in) :: x
-         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
-
-c     outputs
-c     -----------
-
-         DOUBLE PRECISION,INTENT(out) :: Kthmix
-
-c     local
-c    ------------
-
-         DOUBLE PRECISION x1,x2
-
-         DOUBLE PRECISION  Tc1, Tc2, Pc1, Pc2
-
-         DOUBLE PRECISION  A12, A11, A22, A21
-
-         DOUBLE PRECISION  Gamma1, Gamma2, M1, M2
-         DOUBLE PRECISION  lambda_trans1, lambda_trans2,epsilon
-
-         DOUBLE PRECISION  kco2, kn2
-
-      x1 = x
-      x2 = 1d0 - x
-
-      M1 = mco2
-      M2 = mn2
-
-      Tc1 =  304.1282 !(Scalabrin et al. 2006)
-      Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)
-
-      Pc1 =  73.773   ! (bars)
-      Pc2 =  33.958   ! (bars)
-    
-      Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
-      Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)
-
-c Translational conductivities
-
-      lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
-     &                                                          /Gamma1
-
-      lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) )
-     &                                                          /Gamma2
-      
-c     Coefficient of Mason and Saxena
-      epsilon = 1.
-
-      A11 = 1.
-	 
-      A22 = 1.
-
-      A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
-     &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
-
-      A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)*
-     &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
-
-c     INDIVIDUAL COND.
-
-         call KthCO2Scalab(kco2,T,rho)
-         call KthN2LemJac(kn2,T,rho)
-
-c     MIXTURE COND.
-        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
-        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
-
-        END
-
-c======================================================================
-         subroutine KthN2LemJac(kthn2,T,rho)
-c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
-cWITH viscosity
-c======================================================================
-
-       implicit none
-
-        include "microphys.h"
-c        include "microphysCO2.h"
-
-
-c      arguments
-c     -----------
-         
-         REAL,INTENT(in) :: T
-         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
-
-c     outputs
-c     -----------
-
-         DOUBLE PRECISION,INTENT(out) :: kthn2
-
-c     local
-c    ------------
-
-        DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
-        DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
-        DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
-        DOUBLE PRECISION d4,d5,d6,d7,d8,d9
-        DOUBLE PRECISION l4,l5,l6,l7,l8,l9
-        DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9
-        DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9
-
-        DOUBLE PRECISION Tc,rhoc
-
-        DOUBLE PRECISION tau, delta
-
-        DOUBLE PRECISION visco
-
-        DOUBLE PRECISION k1, k2
-
-         N1 = 1.511d0
-         N2 = 2.117d0
-         N3 = -3.332d0
-
-         N4 = 8.862
-         N5 = 31.11
-         N6 = -73.13
-         N7 = 20.03
-         N8 = -0.7096
-         N9 = 0.2672
-
-         t2 = -1.0d0
-         t3 = -0.7d0
-         t4 = 0.0d0
-         t5 = 0.03
-         t6 = 0.2
-         t7 = 0.8
-         t8 = 0.6
-         t9 = 1.9
-   
-         d4 =  1.
-         d5 =  2.
-         d6 =  3.
-         d7 =  4.
-         d8 =  8.
-         d9 = 10.
-   
-         l4 = 0. 
-         gamma4 = 0.
-        
-         l5 = 0. 
-         gamma5 = 0.
-        
-         l6 = 1. 
-         gamma6 = 1.
-
-         l7 = 2. 
-         gamma7 = 1.
-
-         l8 = 2. 
-         gamma8 = 1.
-
-         l9 = 2. 
-         gamma9 = 1.
-
-c----------------------------------------------------------------------           
-         call viscoN2(T,visco)  !! v given in microPa.s
-      
-         Tc   = 126.192d0
-         rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3
-
-         tau  = Tc / T
-         delta = rho/rhoc 
-
-         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1     
-c--------- residual thermal conductivity
-
-         k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
-     &  +     N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) 
-     &  +     N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) 
-     &  +     N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) 
-     &  +     N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) 
-     &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) 
-
-         kthn2 = k1 + k2
-
-         END
-
-
-c======================================================================
-
-         subroutine viscoN2(T,visco)
-
-c        Compute viscosity of N2 (Lemmon and Jacobsen, 2003)
-
-c======================================================================
-
-         implicit none
-
-       include "microphys.h"
-c       include "microphysCO2.h"
-c      arguments
-c     -----------
-         
-      REAL,INTENT(in) :: T
-
-c     outputs
-c     -----------
-
-      DOUBLE PRECISION,INTENT(out) :: visco
-
-
-c     local
-c    ------------
-
-      DOUBLE PRECISION a0,a1,a2,a3,a4 
-      DOUBLE PRECISION Tstar,factor,sigma,M2
-      DOUBLE PRECISION RGCS
-      
-
-c----------------------------------------------------------------------  
-
-  
-      factor = 98.94   ! (K)   
-  
-      sigma  = 0.3656  ! (nm)
-  
-      a0 =  0.431
-      a1 = -0.4623
-      a2 =  0.08406
-      a3 =  0.005341
-      a4 = -0.00331
-  
-      M2 = mn2 * 1.00e3   !!! to g.mol-1
-  
-      Tstar = T*1./factor
-
-      RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + 
-     &                a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. )
-  
-  
-      visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS )  !!! microPa.s
-
-
-      RETURN
-
-      END
-
-
-c======================================================================
-
-         subroutine KthCO2Scalab(kthco2,T,rho)
-
-c        Compute thermal cond of CO2 (Scalabrin et al. 2006)
-
-c======================================================================
-
-         implicit none
-
-
-
-c      arguments
-c     -----------
-
-      REAL,INTENT(in) :: T
-      DOUBLE PRECISION,INTENT(in) :: rho
-
-c     outputs
-c     -----------
-
-      DOUBLE PRECISION,INTENT(out) :: kthco2
-
-c     LOCAL
-c     -----------
-
-      DOUBLE PRECISION Tc,Pc,rhoc, Lambdac
-      
-      DOUBLE PRECISION Tr, rhor, k1, k2
-
-      DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
-      DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
-      DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
-
-      Tc   = 304.1282   !(K) 
-      Pc   = 7.3773e6   !(MPa)
-      rhoc = 467.6      !(kg.m-3)
-      Lambdac = 4.81384 !(mW.m-1K-1)
-  
-      g1 = 0.
-      g2 = 0.
-      g3 = 1.5
-      g4 = 0.0
-      g5 = 1.0
-      g6 = 1.5
-      g7 = 1.5
-      g8 = 1.5
-      g9 = 3.5
-      g10 = 5.5
-
-      h1 = 1.
-      h2 = 5.
-      h3 = 1.
-      h4 = 1.
-      h5 = 2.
-      h6 = 0.
-      h7 = 5.0
-      h8 = 9.0
-      h9 = 0.
-      h10 = 0.
-
-      n1 = 7.69857587
-      n2 = 0.159885811
-      n3 = 1.56918621
-      n4 = -6.73400790
-      n5 = 16.3890156
-      n6 = 3.69415242
-      n7 = 22.3205514
-      n8 = 66.1420950
-      n9 = -0.171779133
-      n10 = 0.00433043347
-
-      Tr   = T/Tc
-      rhor = rho/rhoc
-
-      k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) 
-     &     + n3*Tr**(g3)*rhor**(h3)
-
-      k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5)  
-     &    + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) 
-     &    + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) 
-     &    + n10*Tr**(g10)*rhor**(h10)
-    
-      k2  = exp(-5.*rhor**(2.)) * k2
-                
-      kthco2 = (k1 + k2) *  Lambdac   ! mW
-
-      END
Index: trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F	(revision 2156)
+++ trunk/LMDZ.MARS/libf/phymars/massflowrateco2.F	(revision 2156)
@@ -0,0 +1,540 @@
+c=======================================================================
+      subroutine massflowrateco2(P,T,Sat,Radius,Matm,Ic)
+c
+c     Determination of the mass transfer rate for CO2 condensation & 
+c     sublimation
+c
+c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
+c           particle radius (Radius), molecular mass of the atm (Matm)
+c   output:  MASS FLUX Ic
+c
+c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
+c  
+c
+c   Updates:
+c   --------
+c   December 2017 - C. Listowski - Simplification of the derivation of
+c   massflowrate by using explicit formula for surface temperature,
+c   No Newton-Raphson routine anymore- see comment at relevant line
+c=======================================================================
+      USE comcstfi_h, ONLY: pi
+
+      implicit none
+
+      include "microphys.h"
+
+c   arguments: INPUT
+c   ----------
+      REAL,INTENT(in) :: T,Matm
+      REAL*8,INTENT(in) :: SAT
+      REAL,INTENT(in) :: P
+      DOUBLE PRECISION,INTENT(in) :: Radius
+c   arguments: OUTPUT
+c   ----------
+      DOUBLE PRECISION,INTENT(out) ::   Ic
+c   Local Variables
+c   ----------
+      DOUBLE PRECISION   Tsurf
+      DOUBLE PRECISION   C0,C1,C2
+      DOUBLE PRECISION   kmix,Lsub,cond
+      DOUBLE PRECISION   Ak
+
+      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2,Ak)
+   
+      Tsurf = 1./C1*dlog(Sat/Ak) + T 
+
+      !Note by CL - dec 2017 (see also technical note)
+      !The above is a simplified version of Tsurf
+      !compared to the one used by Listowski et al. 2014 (Ta), where a
+      !Newton-Raphson routine must be used. Approximations made by
+      !considering the orders of magnitude of the different factors lead to
+      !simplification of the equation 5 of Listowski et al. (2014).
+      !The error compared to the exact value determined by NR iterations
+      !is less than 0.6% for all sizes, pressures, supersaturations
+      !relevant to present Mars. Should also be ok for most conditions
+      !in ancient Mars (However, needs to be double-cheked, as explained in
+      !(Listowski et al. 2013,JGR)
+
+      cond = 4.*pi*Radius*kmix
+      Ic = cond*(Tsurf-T)/Lsub
+  
+      END
+
+c********************************************************************************
+
+      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2,Ak)
+
+c********************************************************************************
+c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
+      use tracer_mod, only: rho_ice_co2
+      USE comcstfi_h, ONLY: pi
+
+      implicit none
+      include "microphys.h"
+
+c   arguments: INPUT
+c   ----------------
+      REAL,INTENT(in) :: P
+      REAL,INTENT(in) :: T
+      REAL*8,INTENT(in) :: S
+      DOUBLE PRECISION,INTENT(in) :: rc
+      REAL,INTENT(in) :: Matm !g.mol-1 ( = mmean(ig,l) )
+c   arguments: OUTPUT
+c   ----------
+      DOUBLE PRECISION,INTENT(out) ::  C0,C1,C2
+      DOUBLE PRECISION,INTENT(out) ::  kmix,Lsub
+
+c   local:
+c   ------
+      DOUBLE PRECISION Cpatm,Cpn2,Cpco2
+      DOUBLE PRECISION psat, xinf, pco2
+      DOUBLE PRECISION Dv     
+      DOUBLE PRECISION l0,l1,l2,l3,l4            
+      DOUBLE PRECISION knudsen, a, lambda      ! F and S correction
+      DOUBLE PRECISION Ak                       ! kelvin factor    
+      DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th
+
+
+c     DEFINE heat cap. J.kg-1.K-1 and To
+
+      data Cpco2/0.7e3/
+      data Cpn2/1e3/
+
+      kmix = 0d0
+      Lsub = 0d0
+
+      C0 = 0d0
+      C1 = 0d0
+      C2 = 0d0
+
+c     Equilibirum pressure over a flat surface
+      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
+c     Compute transport coefficient
+      pco2 = psat * dble(S)
+c     Latent heat of sublimation if CO2  co2 (J.kg-1)
+c     version Azreg_Ainou (J/kg) :
+      l0=595594.      
+      l1=903.111     
+      l2=-11.5959    
+      l3=0.0528288 
+      l4=-0.000103183
+      Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * 
+     &     dble(T)**3 + l4 * dble(T)**4 ! J/kg
+c     atmospheric density
+      rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
+      rhoatm = rhoatm * 1.00e-3 !kg.m-3
+      call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
+      call  Diffcoeff(P, T, Dv)  
+
+      Dv = Dv * 1.00e-4         !!! cm2.s-1  to m2.s-1
+
+c     ----- FS correction for Diff
+      vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1
+      knudsen = 3*Dv / (vthco2 * rc) 
+      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
+      Dv      = Dv / (1. + lambda * knudsen)
+c     ----- FS correction for Kth 
+      vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg 
+      Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1
+      lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
+     &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer
+      knudsen = lpmt / rc
+      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
+      kmix    = kmix /  (1. + lambda * knudsen)
+c     --------------------- ASSIGN coeff values for FUNCTION
+      xinf = dble(S) * psat / dble(P)
+      Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
+      C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
+     &     (rgp*dble(T)))
+      C1 = Lsub*mco2/(rgp*dble(T)**2)
+      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
+
+      END
+
+
+c======================================================================
+      subroutine Diffcoeff(P, T, Diff)
+c     Compute diffusion coefficient CO2/N2
+c     cited in Ilona's lecture - from Reid et al. 1987
+c======================================================================
+       IMPLICIT NONE
+
+       include "microphys.h"
+
+c      arguments
+c     -----------
+      
+      REAL,INTENT(in) :: P
+      REAL,INTENT(in) :: T
+      
+c     output
+c     -----------
+
+      DOUBLE PRECISION,INTENT(out) :: Diff
+    
+c      local
+c     -----------
+
+      REAL Pbar                     !!! has to be in bar for the formula
+      DOUBLE PRECISION dva, dvb, Mab  ! Mab has to be in g.mol-1
+        
+        Pbar = P * 1d-5
+    
+        Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000.
+
+  	dva = 26.9        ! diffusion volume of CO2,  Reid et al. 1987 (cited in Ilona's lecture)
+  	dvb = 18.5        ! diffusion volume of N2
+    
+        Diff  = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) 
+     &       * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1
+  
+       RETURN
+
+       END
+
+
+c======================================================================
+
+         subroutine KthMixNEW(Kthmix,T,x,rho)
+
+c        Compute thermal conductivity of CO2/N2 mixture
+c         (***WITHOUT*** USE OF VISCOSITY)
+
+c          (Mason & Saxena, 1958 - Wassiljeva 1904)
+c======================================================================
+
+       implicit none
+       
+       include "microphys.h"
+c      arguments
+c     -----------
+         
+         REAL,INTENT(in) :: T
+         DOUBLE PRECISION,INTENT(in) :: x
+         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
+
+c     outputs
+c     -----------
+
+         DOUBLE PRECISION,INTENT(out) :: Kthmix
+
+c     local
+c    ------------
+
+         DOUBLE PRECISION x1,x2
+
+         DOUBLE PRECISION  Tc1, Tc2, Pc1, Pc2
+
+         DOUBLE PRECISION  A12, A11, A22, A21
+
+         DOUBLE PRECISION  Gamma1, Gamma2, M1, M2
+         DOUBLE PRECISION  lambda_trans1, lambda_trans2,epsilon
+
+         DOUBLE PRECISION  kco2, kn2
+
+      x1 = x
+      x2 = 1d0 - x
+
+      M1 = mco2
+      M2 = mn2
+
+      Tc1 =  304.1282 !(Scalabrin et al. 2006)
+      Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)
+
+      Pc1 =  73.773   ! (bars)
+      Pc2 =  33.958   ! (bars)
+    
+      Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
+      Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)
+
+c Translational conductivities
+
+      lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
+     &                                                          /Gamma1
+
+      lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) )
+     &                                                          /Gamma2
+      
+c     Coefficient of Mason and Saxena
+      epsilon = 1.
+
+      A11 = 1.
+	 
+      A22 = 1.
+
+      A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
+     &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
+
+      A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)*
+     &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
+
+c     INDIVIDUAL COND.
+
+         call KthCO2Scalab(kco2,T,rho)
+         call KthN2LemJac(kn2,T,rho)
+
+c     MIXTURE COND.
+        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
+        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
+
+        END
+
+c======================================================================
+         subroutine KthN2LemJac(kthn2,T,rho)
+c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
+cWITH viscosity
+c======================================================================
+
+       implicit none
+
+        include "microphys.h"
+c        include "microphysCO2.h"
+
+
+c      arguments
+c     -----------
+         
+         REAL,INTENT(in) :: T
+         DOUBLE PRECISION,INTENT(in) :: rho !kg.m-3
+
+c     outputs
+c     -----------
+
+         DOUBLE PRECISION,INTENT(out) :: kthn2
+
+c     local
+c    ------------
+
+        DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
+        DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
+        DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
+        DOUBLE PRECISION d4,d5,d6,d7,d8,d9
+        DOUBLE PRECISION l4,l5,l6,l7,l8,l9
+        DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9
+        DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9
+
+        DOUBLE PRECISION Tc,rhoc
+
+        DOUBLE PRECISION tau, delta
+
+        DOUBLE PRECISION visco
+
+        DOUBLE PRECISION k1, k2
+
+         N1 = 1.511d0
+         N2 = 2.117d0
+         N3 = -3.332d0
+
+         N4 = 8.862
+         N5 = 31.11
+         N6 = -73.13
+         N7 = 20.03
+         N8 = -0.7096
+         N9 = 0.2672
+
+         t2 = -1.0d0
+         t3 = -0.7d0
+         t4 = 0.0d0
+         t5 = 0.03
+         t6 = 0.2
+         t7 = 0.8
+         t8 = 0.6
+         t9 = 1.9
+   
+         d4 =  1.
+         d5 =  2.
+         d6 =  3.
+         d7 =  4.
+         d8 =  8.
+         d9 = 10.
+   
+         l4 = 0. 
+         gamma4 = 0.
+        
+         l5 = 0. 
+         gamma5 = 0.
+        
+         l6 = 1. 
+         gamma6 = 1.
+
+         l7 = 2. 
+         gamma7 = 1.
+
+         l8 = 2. 
+         gamma8 = 1.
+
+         l9 = 2. 
+         gamma9 = 1.
+
+c----------------------------------------------------------------------           
+         call viscoN2(T,visco)  !! v given in microPa.s
+      
+         Tc   = 126.192d0
+         rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3
+
+         tau  = Tc / T
+         delta = rho/rhoc 
+
+         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1     
+c--------- residual thermal conductivity
+
+         k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
+     &  +     N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) 
+     &  +     N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) 
+     &  +     N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) 
+     &  +     N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) 
+     &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) 
+
+         kthn2 = k1 + k2
+
+         END
+
+
+c======================================================================
+
+         subroutine viscoN2(T,visco)
+
+c        Compute viscosity of N2 (Lemmon and Jacobsen, 2003)
+
+c======================================================================
+
+         implicit none
+
+       include "microphys.h"
+c       include "microphysCO2.h"
+c      arguments
+c     -----------
+         
+      REAL,INTENT(in) :: T
+
+c     outputs
+c     -----------
+
+      DOUBLE PRECISION,INTENT(out) :: visco
+
+
+c     local
+c    ------------
+
+      DOUBLE PRECISION a0,a1,a2,a3,a4 
+      DOUBLE PRECISION Tstar,factor,sigma,M2
+      DOUBLE PRECISION RGCS
+      
+
+c----------------------------------------------------------------------  
+
+  
+      factor = 98.94   ! (K)   
+  
+      sigma  = 0.3656  ! (nm)
+  
+      a0 =  0.431
+      a1 = -0.4623
+      a2 =  0.08406
+      a3 =  0.005341
+      a4 = -0.00331
+  
+      M2 = mn2 * 1.00e3   !!! to g.mol-1
+  
+      Tstar = T*1./factor
+
+      RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + 
+     &                a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. )
+  
+  
+      visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS )  !!! microPa.s
+
+
+      RETURN
+
+      END
+
+
+c======================================================================
+
+         subroutine KthCO2Scalab(kthco2,T,rho)
+
+c        Compute thermal cond of CO2 (Scalabrin et al. 2006)
+
+c======================================================================
+
+         implicit none
+
+
+
+c      arguments
+c     -----------
+
+      REAL,INTENT(in) :: T
+      DOUBLE PRECISION,INTENT(in) :: rho
+
+c     outputs
+c     -----------
+
+      DOUBLE PRECISION,INTENT(out) :: kthco2
+
+c     LOCAL
+c     -----------
+
+      DOUBLE PRECISION Tc,Pc,rhoc, Lambdac
+      
+      DOUBLE PRECISION Tr, rhor, k1, k2
+
+      DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
+      DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
+      DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
+
+      Tc   = 304.1282   !(K) 
+      Pc   = 7.3773e6   !(MPa)
+      rhoc = 467.6      !(kg.m-3)
+      Lambdac = 4.81384 !(mW.m-1K-1)
+  
+      g1 = 0.
+      g2 = 0.
+      g3 = 1.5
+      g4 = 0.0
+      g5 = 1.0
+      g6 = 1.5
+      g7 = 1.5
+      g8 = 1.5
+      g9 = 3.5
+      g10 = 5.5
+
+      h1 = 1.
+      h2 = 5.
+      h3 = 1.
+      h4 = 1.
+      h5 = 2.
+      h6 = 0.
+      h7 = 5.0
+      h8 = 9.0
+      h9 = 0.
+      h10 = 0.
+
+      n1 = 7.69857587
+      n2 = 0.159885811
+      n3 = 1.56918621
+      n4 = -6.73400790
+      n5 = 16.3890156
+      n6 = 3.69415242
+      n7 = 22.3205514
+      n8 = 66.1420950
+      n9 = -0.171779133
+      n10 = 0.00433043347
+
+      Tr   = T/Tc
+      rhor = rho/rhoc
+
+      k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) 
+     &     + n3*Tr**(g3)*rhor**(h3)
+
+      k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5)  
+     &    + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) 
+     &    + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) 
+     &    + n10*Tr**(g10)*rhor**(h10)
+    
+      k2  = exp(-5.*rhor**(2.)) * k2
+                
+      kthco2 = (k1 + k2) *  Lambdac   ! mW
+
+      END
Index: trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F	(revision 2155)
+++ 	(revision )
@@ -1,177 +1,0 @@
-*******************************************************
-*                                                     *
-      subroutine nucleaCO2(pco2,temp,sat,n_ccn,nucrate,
-     &           n_ccn_h2oice,rad_h2oice,nucrate_h2oice,
-     &           vo2co2)
-      USE comcstfi_h
-
-      implicit none
-*                                                     *
-*   This subroutine computes the nucleation rate      *
-*   as given in Pruppacher & Klett (1978) in the      *
-*   case of water ice forming on a solid substrate.   *
-*     Definition refined by Keese (jgr,1989)          *
-*   Authors: F. Montmessin                            *
-*     Adapted for the LMD/GCM by J.-B. Madeleine      *
-*     (October 2011)                                  *
-*     Optimisation by A. Spiga (February 2012)        * 
-*   CO2 nucleation routine dev. by Constantino        *
-*     Listowski and Joachim Audouard (2016-2017),     *
-*     adapted from the water ice nucleation     
-* It computes two different nucleation rates : one 
-* on the dust CCN distribution and the other one on
-* the water ice particles distribution
-*******************************************************
- ! nucrate = output
- ! nucrate_h2o en sortie aussi : 
-!nucleation sur dust et h2o separement ici
-
-      include "microphys.h"
-      include "callkeys.h"
-
-c     Inputs
-      DOUBLE PRECISION pco2,sat,vo2co2
-      DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld)
-      REAL temp !temperature
-c     Output
-      DOUBLE PRECISION nucrate(nbinco2_cld)
-      DOUBLE PRECISION nucrate_h2oice(nbinco2_cld) ! h2o as substrate
-      double precision rad_h2oice(nbinco2_cld) ! h2o ice grid (as substrate)
-
-c     Local variables
-      DOUBLE PRECISION nco2
-      DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
-      DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
-      DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
-      DOUBLE PRECISION fshapeco2   ! function defined at the end of the file
-      DOUBLE PRECISION deltaf      
-      double precision mtetalocal,mtetalocalh ! local mteta in double precision
-      double precision fshapeco2simple,zefshapeco2
-      integer i
-c     *************************************************
-
-      mtetalocal = dble(mtetaco2)  !! use mtetalocal for better performance
-      mtetalocalh=dble(mteta)
-
-
-      IF (sat .gt. 1.) THEN    ! minimum condition to activate nucleation
-
-        nco2   = pco2 / kbz / temp
-        rstar  = 2. * sigco2 * vo2co2 / (kbz*temp*dlog(sat))
-        gstar  = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2)
-        
-       fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
-     &                   / 4.
-
-c       Loop over size bins
-        do i=1,nbinco2_cld
-c            write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i),
-c     &          n_ccn(i)
-
-          if ( n_ccn(i) .lt. 1e-10 ) then
-c           no dust, no need to compute nucleation!
-            nucrate(i)=0.
-c            goto 210
-c          endif
-          else
-            if (rad_cldco2(i).gt.3000.*rstar) then
-              zefshapeco2 = fshapeco2simple
-            else
-             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
-            endif
-
-            fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
-     &             zefshapeco2
-            deltaf = (2.*desorpco2-surfdifco2-fistar)/
-     &             (kbz*temp)
-            deltaf = min( max(deltaf, -100.d0), 100.d0)
-
-            if (deltaf.eq.-100.) then
-                nucrate(i) = 0.
-            else
-                nucrate(i)= dble(sqrt ( fistar /
-     &               (3.*pi*kbz*temp*(gstar*gstar)) )
-     &                  * kbz * temp * rstar
-     &                  * rstar * 4. * pi
-     &                  * ( nco2*rad_cldco2(i) )
-     &                  * ( nco2*rad_cldco2(i) )
-     &                  / ( zefshapeco2 * nusco2 * m0co2 )
-     &                  * dexp (deltaf))
-
-            
-            endif
-          endif ! if n_ccn(i) .lt. 1e-10
-
-          if (co2useh2o) then
-
-            if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
-c           no H2O ice, no need to compute nucleation!
-               nucrate_h2oice(i)=0.
-            else   
-               if (rad_h2oice(i).gt.3000.*rstar) then
-                 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
-     &                (1.-mtetalocalh) / 4.
-               else  ! same m for dust/h2o ice
-                 zefshapeco2 = fshapeco2(mtetalocalh,
-     &                            (rad_h2oice(i)/rstar))
-               endif
-
-               fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
-     &             zefshapeco2
-              deltaf = (2.*desorpco2-surfdifco2-fistar)/
-     &             (kbz*temp)
-              deltaf = min( max(deltaf, -100.d0), 100.d0)
-
-              if (deltaf.eq.-100.) then
-                  nucrate_h2oice(i) = 0.
-              else
-                  nucrate_h2oice(i)= dble(sqrt ( fistar /
-     &               (3.*pi*kbz*temp*(gstar*gstar)) )
-     &                  * kbz * temp * rstar
-     &                  * rstar * 4. * pi
-     &                  * ( nco2*rad_h2oice(i) )
-     &                  * ( nco2*rad_h2oice(i) )
-     &                  / ( zefshapeco2 * nusco2 * m0co2 )
-     &                  * dexp (deltaf))
-              endif
-            endif
-          endif
-        enddo
-
-      ELSE ! parcelle d'air non saturée
-
-        do i=1,nbinco2_cld
-          nucrate(i) = 0.
-          nucrate_h2oice(i) = 0.
-        enddo
-
-      ENDIF ! if (sat .gt. 1.) 
-
-      end
-
-*********************************************************
-      double precision function fshapeco2(cost,rap)
-      implicit none
-*        function computing the f(m,x) factor           *
-* related to energy required to form a critical embryo  *
-*********************************************************
-
-      double precision cost,rap
-      double precision yeah
-
-          !! PHI
-          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
-          !! FSHAPECO2 = TERM A
-          fshapeco2 = (1.-cost*rap) / yeah
-          fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2
-          fshapeco2 = 1. + fshapeco2
-          !! ... + TERM B
-          yeah = (rap-cost)/yeah
-          fshapeco2 = fshapeco2 + 
-     & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
-          !! ... + TERM C 
-          fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.)
-          !! FACTOR 1/2
-          fshapeco2 = 0.5*fshapeco2
-
-      end
Index: trunk/LMDZ.MARS/libf/phymars/nucleaco2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nucleaco2.F	(revision 2156)
+++ trunk/LMDZ.MARS/libf/phymars/nucleaco2.F	(revision 2156)
@@ -0,0 +1,177 @@
+*******************************************************
+*                                                     *
+      subroutine nucleaco2(pco2,temp,sat,n_ccn,nucrate,
+     &           n_ccn_h2oice,rad_h2oice,nucrate_h2oice,
+     &           vo2co2)
+      USE comcstfi_h
+
+      implicit none
+*                                                     *
+*   This subroutine computes the nucleation rate      *
+*   as given in Pruppacher & Klett (1978) in the      *
+*   case of water ice forming on a solid substrate.   *
+*     Definition refined by Keese (jgr,1989)          *
+*   Authors: F. Montmessin                            *
+*     Adapted for the LMD/GCM by J.-B. Madeleine      *
+*     (October 2011)                                  *
+*     Optimisation by A. Spiga (February 2012)        * 
+*   CO2 nucleation routine dev. by Constantino        *
+*     Listowski and Joachim Audouard (2016-2017),     *
+*     adapted from the water ice nucleation     
+* It computes two different nucleation rates : one 
+* on the dust CCN distribution and the other one on
+* the water ice particles distribution
+*******************************************************
+ ! nucrate = output
+ ! nucrate_h2o en sortie aussi : 
+!nucleation sur dust et h2o separement ici
+
+      include "microphys.h"
+      include "callkeys.h"
+
+c     Inputs
+      DOUBLE PRECISION pco2,sat,vo2co2
+      DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld)
+      REAL temp !temperature
+c     Output
+      DOUBLE PRECISION nucrate(nbinco2_cld)
+      DOUBLE PRECISION nucrate_h2oice(nbinco2_cld) ! h2o as substrate
+      double precision rad_h2oice(nbinco2_cld) ! h2o ice grid (as substrate)
+
+c     Local variables
+      DOUBLE PRECISION nco2
+      DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
+      DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
+      DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
+      DOUBLE PRECISION fshapeco2   ! function defined at the end of the file
+      DOUBLE PRECISION deltaf      
+      double precision mtetalocal,mtetalocalh ! local mteta in double precision
+      double precision fshapeco2simple,zefshapeco2
+      integer i
+c     *************************************************
+
+      mtetalocal = dble(mtetaco2)  !! use mtetalocal for better performance
+      mtetalocalh=dble(mteta)
+
+
+      IF (sat .gt. 1.) THEN    ! minimum condition to activate nucleation
+
+        nco2   = pco2 / kbz / temp
+        rstar  = 2. * sigco2 * vo2co2 / (kbz*temp*dlog(sat))
+        gstar  = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2)
+        
+       fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
+     &                   / 4.
+
+c       Loop over size bins
+        do i=1,nbinco2_cld
+c            write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i),
+c     &          n_ccn(i)
+
+          if ( n_ccn(i) .lt. 1e-10 ) then
+c           no dust, no need to compute nucleation!
+            nucrate(i)=0.
+c            goto 210
+c          endif
+          else
+            if (rad_cldco2(i).gt.3000.*rstar) then
+              zefshapeco2 = fshapeco2simple
+            else
+             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
+            endif
+
+            fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
+     &             zefshapeco2
+            deltaf = (2.*desorpco2-surfdifco2-fistar)/
+     &             (kbz*temp)
+            deltaf = min( max(deltaf, -100.d0), 100.d0)
+
+            if (deltaf.eq.-100.) then
+                nucrate(i) = 0.
+            else
+                nucrate(i)= dble(sqrt ( fistar /
+     &               (3.*pi*kbz*temp*(gstar*gstar)) )
+     &                  * kbz * temp * rstar
+     &                  * rstar * 4. * pi
+     &                  * ( nco2*rad_cldco2(i) )
+     &                  * ( nco2*rad_cldco2(i) )
+     &                  / ( zefshapeco2 * nusco2 * m0co2 )
+     &                  * dexp (deltaf))
+
+            
+            endif
+          endif ! if n_ccn(i) .lt. 1e-10
+
+          if (co2useh2o) then
+
+            if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
+c           no H2O ice, no need to compute nucleation!
+               nucrate_h2oice(i)=0.
+            else   
+               if (rad_h2oice(i).gt.3000.*rstar) then
+                 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
+     &                (1.-mtetalocalh) / 4.
+               else  ! same m for dust/h2o ice
+                 zefshapeco2 = fshapeco2(mtetalocalh,
+     &                            (rad_h2oice(i)/rstar))
+               endif
+
+               fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * 
+     &             zefshapeco2
+              deltaf = (2.*desorpco2-surfdifco2-fistar)/
+     &             (kbz*temp)
+              deltaf = min( max(deltaf, -100.d0), 100.d0)
+
+              if (deltaf.eq.-100.) then
+                  nucrate_h2oice(i) = 0.
+              else
+                  nucrate_h2oice(i)= dble(sqrt ( fistar /
+     &               (3.*pi*kbz*temp*(gstar*gstar)) )
+     &                  * kbz * temp * rstar
+     &                  * rstar * 4. * pi
+     &                  * ( nco2*rad_h2oice(i) )
+     &                  * ( nco2*rad_h2oice(i) )
+     &                  / ( zefshapeco2 * nusco2 * m0co2 )
+     &                  * dexp (deltaf))
+              endif
+            endif
+          endif
+        enddo
+
+      ELSE ! parcelle d'air non saturée
+
+        do i=1,nbinco2_cld
+          nucrate(i) = 0.
+          nucrate_h2oice(i) = 0.
+        enddo
+
+      ENDIF ! if (sat .gt. 1.) 
+
+      end
+
+*********************************************************
+      double precision function fshapeco2(cost,rap)
+      implicit none
+*        function computing the f(m,x) factor           *
+* related to energy required to form a critical embryo  *
+*********************************************************
+
+      double precision cost,rap
+      double precision yeah
+
+          !! PHI
+          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
+          !! FSHAPECO2 = TERM A
+          fshapeco2 = (1.-cost*rap) / yeah
+          fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2
+          fshapeco2 = 1. + fshapeco2
+          !! ... + TERM B
+          yeah = (rap-cost)/yeah
+          fshapeco2 = fshapeco2 + 
+     & rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
+          !! ... + TERM C 
+          fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.)
+          !! FACTOR 1/2
+          fshapeco2 = 0.5*fshapeco2
+
+      end
