Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 2349)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 2349)
@@ -0,0 +1,1022 @@
+      MODULE co2cloud_mod
+
+      IMPLICIT NONE 
+
+      DOUBLE PRECISION,allocatable,save :: mem_Mccn_co2(:,:) ! Memory of CCN mass of H2O and dust used by CO2
+      DOUBLE PRECISION,allocatable,save :: mem_Mh2o_co2(:,:) ! Memory of H2O mass integred into CO2 crystal
+      DOUBLE PRECISION,allocatable,save :: mem_Nccn_co2(:,:) ! Memory of CCN number of H2O and dust used by CO2
+
+      CONTAINS
+
+      SUBROUTINE co2cloud(ngrid,nlay,ptimestep, 
+     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
+     &                pq,pdq,pdqcloudco2,pdtcloudco2,
+     &                nq,tau,tauscaling,rdust,rice,riceco2,nuice,
+     &                rsedcloudco2,rhocloudco2,
+     &                rsedcloud,rhocloud,pzlev,pdqs_sedco2,
+     &                pdu,pu,pcondicea)
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+      use dimradmars_mod, only: naerkind
+      USE comcstfi_h, only: pi, g, cpp
+      USE updaterad, only: updaterice_microco2, updaterice_micro,
+     &                     updaterdust
+      use conc_mod, only: mmean,rnew
+      use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
+     &     igcm_dust_mass, igcm_dust_number,igcm_h2o_ice,
+     &     igcm_ccn_mass,igcm_ccn_number,
+     &     igcm_ccnco2_mass, igcm_ccnco2_number,
+     &     rho_dust, nuiceco2_sed, nuiceco2_ref,
+     &     rho_ice_co2,r3n_q,rho_ice,nuice_sed
+      USE newsedim_mod, ONLY: newsedim
+      USE datafile_mod, ONLY: datadir
+      USE improvedco2clouds_mod, ONLY: improvedco2clouds
+
+      IMPLICIT NONE
+
+      include "callkeys.h"
+      include "microphys.h"
+
+c=======================================================================
+c CO2 clouds formation
+c
+c  There is a time loop specific to cloud formation 
+c  due to timescales smaller than the GCM integration timestep.
+c  microphysics subroutine is improvedCO2clouds.F
+c  the microphysics time step is a fraction of the physical one
+c  the integer imicroco2 must be set in callphys.def  
+c
+c  The co2 clouds tracers (co2_ice, ccn mass and concentration) are 
+c  sedimented at each microtimestep. pdqs_sedco2 keeps track of the 
+c  CO2 flux at the surface 
+c
+c  Authors: 09/2016 Joachim Audouard & Constantino Listowski 
+c  Adaptation of the water ice clouds scheme (with specific microphysics)
+c  of Montmessin, Navarro & al.
+c
+c 07/2017 J.Audouard
+c Several logicals and integer must be set to .true. in callphys.def 
+c if not, default values are .false.
+c co2clouds=.true. call this routine
+c co2useh2o=.true. allow the use of water ice particles as CCN for CO2
+c meteo_flux=.true. supply meteoritic particles 
+c CLFvaryingCO2=.true. allows a subgrid temperature distribution
+c of amplitude spantCO2(=integer in callphys.def, typically 3)
+c satindexco2=.true. allows the filtering out of the sub-grid T distribution
+c                    if the GW saturates in the column. Based on Spiga et al
+c                    2012
+c                    An index is computed for the column, and the sub-grid T 
+c                    distribution is applied if the index remains < 0.1
+c                    setting to .false. applies the sub-grid T everywhere.
+c                    default value is .true., only applies if 
+c                    CLFvaryingCO2=.true. anyway.
+c imicroco2=50 
+c
+c The subgrid Temperature distribution is modulated (0 or 1) by Spiga et
+c al. (GRL 2012) Saturation Index to account for GW propagation or 
+c dissipation upwards.
+c
+c 4D and column opacities are computed using Qext values at 1µm. 
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   arguments:
+c   -------------
+
+      INTEGER, INTENT(IN) :: ngrid,nlay
+      REAL, INTENT(IN) :: ptimestep            ! pas de temps physique (s)
+      REAL, INTENT(IN) :: pplev(ngrid,nlay+1)   ! Inter-layer pressures (Pa)
+      REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! mid-layer pressures (Pa)
+      REAL, INTENT(IN) :: pdpsrf(ngrid)         ! tendency on surface pressure
+      REAL, INTENT(IN) :: pzlay(ngrid,nlay)     ! altitude at the middle of the layers
+      REAL, INTENT(IN) :: pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
+      REAL, INTENT(IN) :: pdt(ngrid,nlay)       ! tendency on temperature from other parametrizations
+      real, INTENT(IN) :: pq(ngrid,nlay,nq)     ! tracers (kg/kg)
+      real, INTENT(IN) :: pdq(ngrid,nlay,nq)    ! tendencies before condensation  (kg/kg.s-1)
+      real, intent(OUT) :: pdqcloudco2(ngrid,nlay,nq) ! tendency due to CO2 condensation (kg/kg.s-1)
+      real, intent(OUT) :: pcondicea(ngrid,nlay)
+      real, intent(OUT) :: pdtcloudco2(ngrid,nlay)    ! tendency on temperature due to latent heat
+      INTEGER, INTENT(IN) :: nq                 ! number of tracers
+      REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point
+      REAL, INTENT(IN) :: tauscaling(ngrid)   ! Convertion factor for dust amount
+      REAL, INTENT(OUT) :: rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
+      real, intent(OUT) :: rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
+                                ! used for nucleation of CO2 on ice-coated ccns
+      DOUBLE PRECISION, INTENT(out) :: riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
+                               ! (r_c in montmessin_2004)
+      REAL, INTENT(IN) :: nuice(ngrid,nlay)   ! Estimated effective variance
+                               !   of the size distribution
+      real, intent(OUT) :: rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
+      real, intent(OUT) :: rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
+      real, intent(OUT) :: rsedcloud(ngrid,nlay) ! Water Cloud sedimentation radius
+      real, intent(OUT) :: rhocloud(ngrid,nlay)  ! Water Cloud density (kg.m-3)
+      real, intent(IN) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
+      real, intent(OUT) :: pdqs_sedco2(ngrid) ! CO2 flux at the surface
+      REAL, INTENT(IN) :: pdu(ngrid,nlay),pu(ngrid,nlay) !Zonal Wind: zu=pu+pdu*ptimestep
+
+c   local:
+c   ------
+        
+      ! for time loop
+      INTEGER microstep  ! time subsampling step variable
+      INTEGER, SAVE :: imicroco2     ! time subsampling for coupled water microphysics & sedimentation
+      REAL, SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation
+      
+      ! tendency given by clouds (inside the micro loop)
+      REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
+      REAL subpdtcloudco2(ngrid,nlay)    ! cf. pdtcloud
+
+      ! global tendency (clouds+physics)
+      REAL sum_subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
+      REAL sum_subpdt(ngrid,nlay)         ! cf. pdtcloud
+      real wq(ngrid,nlay+1)  !  ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
+
+      REAL satuco2(ngrid,nlay)  ! co2 satu ratio for output
+      REAL zqsatco2(ngrid,nlay) ! saturation co2
+
+      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
+      DOUBLE PRECISION :: myT   ! temperature scalar for co2 density computation
+
+      INTEGER iq,ig,l,i
+      LOGICAL,SAVE :: firstcall=.true.
+      DOUBLE PRECISION Nccnco2, Niceco2,Nco2,Qccnco2
+      real :: beta ! for sedimentation
+
+      real epaisseur (ngrid,nlay) ! Layer thickness (m)
+      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
+      real ztsed(ngrid,nlay) ! tracers with real-time value in microtimeloop
+      real zqsed(ngrid,nlay,nq)
+      real zqsed0(ngrid,nlay,nq) !For sedimentation tendancy
+      real subpdqsed(ngrid,nlay,nq)
+      real sum_subpdqs_sedco2(ngrid) ! CO2 flux at the surface
+     
+!     What we need for Qext reading and tau computation : size distribution
+      DOUBLE PRECISION vrat_cld ! Volume ratio
+      DOUBLE PRECISION, SAVE :: rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
+      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 boundary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
+      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)
+      REAL, SAVE :: sigma_iceco2 ! Variance of the ice and CCN distributions
+      logical :: file_ok !Qext file reading
+      double precision :: radv(10000),Qextv1mic(10000)
+      double precision, save :: Qext1bins(100)
+      double precision :: Qtemp
+      double precision :: ltemp1(10000),ltemp2(10000)
+      integer :: nelem,lebon1,lebon2
+      integer,parameter :: uQext=555
+      DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
+      DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
+      DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm
+     
+        ! For sub grid T distribution
+
+      REAL zt(ngrid,nlay)       ! local value of temperature
+      REAL :: zq(ngrid, nlay,nq)
+
+      real :: rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
+
+      DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature
+      REAL ::  zqvap(ngrid,nlay)
+      REAL ::  zqice(ngrid,nlay)
+      REAL ::  spant,zdelt ! delta T for the temperature distribution
+      REAL ::  pteff(ngrid, nlay)! effective temperature in the cloud,neb
+      REAL ::  pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
+      REAL ::  co2cloudfrac(ngrid,nlay) ! cloud fraction
+      REAL ::  mincloud ! min cloud frac
+      DOUBLE PRECISION:: rho,zu,NN,gradT !For Saturation Index computation
+      DOUBLE PRECISION :: SatIndex(ngrid,nlay),SatIndexmap(ngrid)
+
+c D. BARDET : sensibility test
+      REAL :: No_dust(ngrid,nlay)
+      REAL :: Mo_dust(ngrid,nlay)
+      
+c      logical :: CLFvaryingCO2
+c     ** un petit test de coherence
+c       --------------------------
+
+      IF (firstcall) THEN
+        if (nq.gt.nqmx) then
+           write(*,*) 'stop in co2cloud (nq.gt.nqmx)!'
+           write(*,*) 'nq=',nq,' nqmx=',nqmx
+           call abort_physic("co2cloud","nq.gt.nqmx",1)
+        endif
+        write(*,*) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2
+        write(*,*) "co2cloud: igcm_co2=",igcm_co2
+        write(*,*) "            igcm_co2_ice=",igcm_co2_ice
+                
+        write(*,*) "time subsampling for microphysic ?"
+#ifdef MESOSCALE
+        imicroco2 = 2
+#else
+        imicroco2 = 30
+#endif
+        call getin_p("imicroco2",imicroco2)
+        write(*,*)"imicroco2 = ",imicroco2
+        
+        microtimestep = ptimestep/real(imicroco2)
+        write(*,*)"Physical timestep is",ptimestep 
+        write(*,*)"CO2 Microphysics timestep is",microtimestep 
+
+c     Compute the size bins of the distribution of CO2 ice particles
+c --> used for opacity calculations
+
+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.
+        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
+c       Volume ratio between two adjacent bins
+   !     vrat_cld 
+        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)
+
+c   read the Qext values 
+        INQUIRE(FILE=TRIM(datadir)//
+     &       '/optprop_co2ice_1mic.dat', EXIST=file_ok)
+        IF (.not. file_ok) THEN 
+           write(*,*) 'file optprop_co2ice_1mic.dat should be in '
+     &          ,trim(datadir)
+           call abort_physic("co2cloud",
+     &          "missing file optprop_co2ice_1mic.dat",1)
+        endif
+!        open(newunit=uQext,file=trim(datadir)//
+        open(unit=uQext,file=trim(datadir)//
+     &       '/optprop_co2ice_1mic.dat'
+     &       ,FORM='formatted')
+        read(uQext,*) !skip 1 line
+        do i=1,10000 
+           read(uQext,'(E11.5)') radv(i)
+        enddo
+        read(uQext,*) !skip 1 line
+        do i=1,10000 
+           read(uQext,'(E11.5)') Qextv1mic(i)
+        enddo
+        close(uQext)
+c   innterpol the Qext values 
+        !rice_out=rad_cldco2
+        do i=1,nbinco2_cld
+           ltemp1=abs(radv(:)-rb_cldco2(i))
+           ltemp2=abs(radv(:)-rb_cldco2(i+1))
+           lebon1=minloc(ltemp1,DIM=1)
+           lebon2=min(minloc(ltemp2,DIM=1),10000)
+           nelem=lebon2-lebon1+1.
+           Qtemp=0d0
+           do l=0,nelem
+              Qtemp=Qtemp+Qextv1mic(min(lebon1+l,10000)) !mean value in the interval
+           enddo
+           Qtemp=Qtemp/nelem
+           Qext1bins(i)=Qtemp
+        enddo
+        Qext1bins(:)=Qext1bins(:)*rad_cldco2(:)*rad_cldco2(:)*pi
+!     The actuall tau computation and output is performed in co2cloud.F
+
+        print*,'--------------------------------------------'
+        print*,'Microphysics co2: size bin-Qext information:'
+        print*,'   i, rad_cldco2(i), Qext1bins(i)'
+        do i=1,nbinco2_cld
+          write(*,'(i3,3x,3(e13.6,4x))') i, rad_cldco2(i),
+     &    Qext1bins(i)
+        enddo
+        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
+        if (CLFvaryingCO2) then
+          write(*,*) 
+          write(*,*) "CLFvaryingCO2 is set to true is callphys.def"
+          write(*,*) "The temperature field is enlarged to +/-",spantCO2 
+          write(*,*) "for the CO2 microphysics "
+        endif
+
+        firstcall=.false.
+      ENDIF                     ! of IF (firstcall)
+c ===========================================================================   
+c     Initialization
+c ===========================================================================
+      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
+      beta=0.85
+      sum_subpdq(1:ngrid,1:nlay,1:nq) = 0
+      sum_subpdt(1:ngrid,1:nlay)      = 0
+      subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
+      subpdtcloudco2(1:ngrid,1:nlay)      = 0
+
+      wq(:,:)=0
+      ! default value if no ice
+      rhocloudco2(1:ngrid,1:nlay) = rho_dust
+      rhocloudco2t(1:ngrid,1:nlay) = rho_dust
+      epaisseur(1:ngrid,1:nlay)=0
+      masse(1:ngrid,1:nlay)=0
+
+      zqsed0(1:ngrid,1:nlay,1:nq)=0
+      sum_subpdqs_sedco2(1:ngrid)=0
+      subpdqsed(1:ngrid,1:nlay,1:nq)=0
+      
+      do  l=1,nlay
+        do ig=1, ngrid
+          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
+       enddo
+      enddo
+
+c ==========================================================================
+c   0.  Representation of sub-grid water ice clouds
+c ==========================================================================
+      IF (CLFvaryingCO2) THEN
+
+         spant=spantCO2         ! delta T for the temprature distribution
+         mincloud=0.1           ! min co2cloudfrac when there is ice  
+         pteff(:,:)=pt(:,:) 
+         co2cloudfrac(:,:)=mincloud
+         
+c Tendencies
+         DO l=1,nlay
+            DO ig=1,ngrid
+               zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
+            ENDDO
+         ENDDO
+         DO l=1,nlay
+            DO ig=1,ngrid
+               DO iq=1,nq
+                  zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep 
+               ENDDO
+            ENDDO
+         ENDDO
+         zqvap=zq(:,:,igcm_co2)
+         zqice=zq(:,:,igcm_co2_ice)
+
+
+         call WRITEDIAGFI(ngrid,"co2cloud_pzlev","pzlev","km",3,
+     &        pzlev)
+         call WRITEDIAGFI(ngrid,"co2cloud_pzlay","pzlay","km",3,
+     &        pzlay)
+         call WRITEDIAGFI(ngrid,"co2cloud_pplay","pplay","Pa",3,
+     &        pplay)
+
+         if (satindexco2) then !logical in callphys.def
+           DO l=12,26
+             ! layers 12 --> 26 ~ 12->85 km
+             DO ig=1,ngrid
+             ! compute N^2 static stability
+             gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
+             NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
+             ! compute absolute value of zonal wind field 
+             zu=abs(pu(ig,l)  + pdu(ig,l)*ptimestep)
+             ! compute background density
+             rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+             !saturation index
+             SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/
+     &                                               (rho*zu*zu*zu))
+             ENDDO
+           ENDDO
+            !Then compute Satindex map 
+            ! layers 12 --> 26 ~ 12->85 km
+           DO ig=1,ngrid
+             SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
+           ENDDO
+
+           call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
+     &         SatIndexmap)
+         else
+           do ig=1,ngrid
+             SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
+           enddo
+         endif ! of if (satindexco2)
+
+!Modulate the DeltaT by GW propagation index : 
+         ! Saturation index S in Spiga 2012 paper 
+         !Assuming like in the paper, 
+         !GW phase speed (stationary waves) c=0 m.s-1
+         !lambdaH =150 km
+         !Fo=7.5e-7 J.m-3
+        
+         CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
+         zdelt=spant  
+         DO ig=1,ngrid
+              
+           IF (SatIndexmap(ig) .le. 0.1) THEN
+             DO l=1,nlay-1
+         
+               IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt) 
+     &             .or. tcond(ig,l) .le. 0 ) THEN !The entire fraction is saturated
+                  pteff(ig,l)=zt(ig,l)
+                  co2cloudfrac(ig,l)=1.
+               ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN ! No saturation at all
+                  pteff(ig,l)=zt(ig,l)-zdelt
+                  co2cloudfrac(ig,l)=mincloud
+               ELSE
+                  co2cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
+     &                 (2.0*zdelt)
+                  pteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Mean temperature of the cloud fraction
+               END IF           !ig if (tcond(ig,l) ...
+               pteff(ig,l)=pteff(ig,l)-pdt(ig,l)*ptimestep
+               IF (co2cloudfrac(ig,l).le. mincloud) THEN
+                  co2cloudfrac(ig,l)=mincloud
+               ELSE IF (co2cloudfrac(ig,l).gt. 1) THEN
+                  co2cloudfrac(ig,l)=1.
+               END IF
+             ENDDO
+           ELSE
+! SatIndex not favorable for GW : leave pt untouched
+             pteff(ig,l)=pt(ig,l)
+             co2cloudfrac(ig,l)=mincloud
+           ENDIF                 ! of if(SatIndexmap...
+         ENDDO ! of DO ig=1,ngrid
+! Totalcloud frac of the column missing here
+c
+c No sub-grid cloud representation (CLFvarying=false)
+      ELSE
+         DO l=1,nlay
+            DO ig=1,ngrid
+               pteff(ig,l)=pt(ig,l)
+            END DO
+         END DO
+      END IF                    ! end if (CLFvaryingco2)
+c =============================================================================
+c microtimestep timeloop for microphysics:
+c 0.Stepped entry for tendancies
+c 1.Compute sedimentation and update tendancies
+c 2.Call co2clouds microphysics
+c 3.Update tendancies 
+c =============================================================================
+      DO microstep=1,imicroco2
+c Temperature tendency subpdt
+        ! If imicro=1 subpdt is the same as pdt
+        DO l=1,nlay
+          DO ig=1,ngrid
+               sum_subpdt(ig,l) = sum_subpdt(ig,l)
+     &              + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
+               sum_subpdq(ig,l,igcm_dust_mass) = 
+     &              sum_subpdq(ig,l,igcm_dust_mass)
+     &              + pdq(ig,l,igcm_dust_mass)
+               sum_subpdq(ig,l,igcm_dust_number) = 
+     &              sum_subpdq(ig,l,igcm_dust_number)
+     &              + pdq(ig,l,igcm_dust_number)
+
+               sum_subpdq(ig,l,igcm_ccnco2_mass) = 
+     &              sum_subpdq(ig,l,igcm_ccnco2_mass)
+     &              + pdq(ig,l,igcm_ccnco2_mass)
+               sum_subpdq(ig,l,igcm_ccnco2_number) = 
+     &              sum_subpdq(ig,l,igcm_ccnco2_number)
+     &              + pdq(ig,l,igcm_ccnco2_number)
+
+               sum_subpdq(ig,l,igcm_co2_ice) = 
+     &              sum_subpdq(ig,l,igcm_co2_ice)
+     &              + pdq(ig,l,igcm_co2_ice)
+               sum_subpdq(ig,l,igcm_co2) = 
+     &              sum_subpdq(ig,l,igcm_co2)
+     &              + pdq(ig,l,igcm_co2)
+c D.BARDET :
+               if (co2useh2o) then
+               sum_subpdq(ig,l,igcm_h2o_ice) = 
+     &              sum_subpdq(ig,l,igcm_h2o_ice)
+     &              + pdq(ig,l,igcm_h2o_ice)
+
+               sum_subpdq(ig,l,igcm_ccn_mass) = 
+     &              sum_subpdq(ig,l,igcm_ccn_mass)
+     &              + pdq(ig,l,igcm_ccn_mass)
+               sum_subpdq(ig,l,igcm_ccn_number) = 
+     &              sum_subpdq(ig,l,igcm_ccn_number)
+     &              + pdq(ig,l,igcm_ccn_number)
+                endif
+          ENDDO
+        ENDDO
+c Effective tracers quantities in the cloud fraction
+        IF (CLFvaryingCO2) THEN     
+            pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
+            pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
+     &           co2cloudfrac(:,:)
+            pqeff(:,:,igcm_ccnco2_number)= 
+     &           pq(:,:,igcm_ccnco2_number)/co2cloudfrac(:,:)
+            pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
+     &           co2cloudfrac(:,:)
+        ELSE
+            pqeff(:,:,:)=pq(:,:,:)
+        END IF  
+       
+c ========================================================================
+c 1.SEDIMENTATION : update tracers, compute parameters,
+c   call to sedimentation routine, update tendancies
+c ========================================================================
+        IF (sedimentation) THEN
+        
+        DO l=1, nlay
+          DO ig=1,ngrid             
+             ztsed(ig,l)=pteff(ig,l)
+     &            +sum_subpdt(ig,l)*microtimestep
+             zqsed(ig,l,:)=pqeff(ig,l,:)
+     &            +sum_subpdq(ig,l,:)*microtimestep
+             rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*
+     &            ztsed(ig,l)-2.87e-6*
+     &            ztsed(ig,l)*ztsed(ig,l))
+             
+             rho_ice_co2=rho_ice_co2T(ig,l)
+             Niceco2=max(zqsed(ig,l,igcm_co2_ice),1.e-30)
+             Nccnco2=max(zqsed(ig,l,igcm_ccnco2_number),
+     &            1.e-30)
+             Qccnco2=max(zqsed(ig,l,igcm_ccnco2_mass),
+     &            1.e-30)
+             call updaterice_microco2(Niceco2,
+     &            Qccnco2,Nccnco2,
+     &            tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
+             if (Niceco2 .le. 1.e-25 
+     &            .or. Nccnco2*tauscaling(ig) .le. 1) THEN
+                riceco2(ig,l)=1.e-9
+             endif
+             rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l)
+     &            ,rho_ice_co2),rho_dust)
+             rsedcloudco2(ig,l)=max(riceco2(ig,l)*
+     &            (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
+     &            riceco2(ig,l))
+          ENDDO
+        ENDDO
+! Gravitational sedimentation       
+        zqsed0(:,:,igcm_co2_ice)=zqsed(:,:,igcm_co2_ice)
+        zqsed0(:,:,igcm_ccnco2_mass)=zqsed(:,:,igcm_ccnco2_mass)
+        zqsed0(:,:,igcm_ccnco2_number)=zqsed(:,:,igcm_ccnco2_number) 
+! We save actualized tracer values to compute sedimentation tendancies
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &     microtimestep,pplev,masse,epaisseur,ztsed,
+     &     rsedcloudco2,rhocloudco2t,
+     &     zqsed(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
+
+! sedim at the surface of co2 ice : keep track of it for physiq_mod
+        do ig=1,ngrid 
+          sum_subpdqs_sedco2(ig)=
+     &         sum_subpdqs_sedco2(ig)+ wq(ig,1)/microtimestep
+        end do
+
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &     microtimestep,pplev,masse,epaisseur,ztsed,
+     &     rsedcloudco2,rhocloudco2t,
+     &     zqsed(:,:,igcm_ccnco2_mass),wq,beta) 
+
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &     microtimestep,pplev,masse,epaisseur,ztsed,
+     &     rsedcloudco2,rhocloudco2t,
+     &     zqsed(:,:,igcm_ccnco2_number),wq,beta) 
+
+        DO l = 1, nlay            !Compute tendencies
+          DO ig=1,ngrid
+            subpdqsed(ig,l,igcm_ccnco2_mass)=
+     &           (zqsed(ig,l,igcm_ccnco2_mass)-
+     &           zqsed0(ig,l,igcm_ccnco2_mass))/microtimestep
+            subpdqsed(ig,l,igcm_ccnco2_number)=
+     &           (zqsed(ig,l,igcm_ccnco2_number)-
+     &           zqsed0(ig,l,igcm_ccnco2_number))/microtimestep
+            subpdqsed(ig,l,igcm_co2_ice)=
+     &           (zqsed(ig,l,igcm_co2_ice)-
+     &           zqsed0(ig,l,igcm_co2_ice))/microtimestep
+          ENDDO
+        ENDDO
+! update subtimestep tendencies with sedimentation input
+        DO l=1,nlay
+         DO ig=1,ngrid
+            sum_subpdq(ig,l,igcm_ccnco2_mass) =
+     &           sum_subpdq(ig,l,igcm_ccnco2_mass)
+     &           +subpdqsed(ig,l,igcm_ccnco2_mass)
+            sum_subpdq(ig,l,igcm_ccnco2_number) =
+     &           sum_subpdq(ig,l,igcm_ccnco2_number)
+     &           +subpdqsed(ig,l,igcm_ccnco2_number)
+            sum_subpdq(ig,l,igcm_co2_ice) =
+     &           sum_subpdq(ig,l,igcm_co2_ice)
+     &           +subpdqsed(ig,l,igcm_co2_ice)
+         ENDDO
+        ENDDO
+        
+        END IF !(end if sedimentation)
+        
+c ==============================================================================
+c      2.  Main call to the cloud schemes:
+c ==============================================================================
+        CALL 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)
+c D. BARDET: sensibility test
+c      call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussiere"
+c     &        ,"part/kg",3,No_dust)
+c      call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussiere"
+c     &        ,"kg/kg ",3,Mo_dust)
+c ==============================================================================
+c      3.  Updating tendencies after cloud scheme:
+c ==============================================================================
+        DO l=1,nlay
+          DO ig=1,ngrid
+               sum_subpdt(ig,l) =
+     &              sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
+
+               sum_subpdq(ig,l,igcm_dust_mass) =
+     &              sum_subpdq(ig,l,igcm_dust_mass)
+     &              + subpdqcloudco2(ig,l,igcm_dust_mass)
+               sum_subpdq(ig,l,igcm_dust_number) =
+     &              sum_subpdq(ig,l,igcm_dust_number)
+     &              + subpdqcloudco2(ig,l,igcm_dust_number)
+
+               sum_subpdq(ig,l,igcm_ccnco2_mass) =
+     &              sum_subpdq(ig,l,igcm_ccnco2_mass)
+     &              + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
+               sum_subpdq(ig,l,igcm_ccnco2_number) =
+     &              sum_subpdq(ig,l,igcm_ccnco2_number)
+     &              + subpdqcloudco2(ig,l,igcm_ccnco2_number)
+
+               sum_subpdq(ig,l,igcm_co2_ice) =
+     &              sum_subpdq(ig,l,igcm_co2_ice)
+     &              + subpdqcloudco2(ig,l,igcm_co2_ice)
+               sum_subpdq(ig,l,igcm_co2) =
+     &              sum_subpdq(ig,l,igcm_co2)
+     &              + subpdqcloudco2(ig,l,igcm_co2)
+c D.BARDET :
+               if (co2useh2o) then
+               sum_subpdq(ig,l,igcm_h2o_ice) =
+     &              sum_subpdq(ig,l,igcm_h2o_ice)
+     &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
+
+               sum_subpdq(ig,l,igcm_ccn_mass) =
+     &              sum_subpdq(ig,l,igcm_ccn_mass)
+     &              + subpdqcloudco2(ig,l,igcm_ccn_mass)
+               sum_subpdq(ig,l,igcm_ccn_number) =
+     &              sum_subpdq(ig,l,igcm_ccn_number)
+     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
+                endif
+          ENDDO
+        ENDDO
+      ENDDO                     ! of DO microstep=1,imicro
+      
+c------------------------------------------------
+c   Compute final tendencies after time loop:
+c------------------------------------------------
+c Condensation/sublimation tendency after clouds scheme (to replace
+c zcondicea in newcondens.F
+
+      DO l=nlay, 1, -1
+        DO ig = 1, ngrid
+           pcondicea(ig,l) = sum_subpdq(ig,l,igcm_co2_ice)
+     &                       /real(imicroco2)
+        ENDDO
+      ENDDO
+
+c CO2 flux at surface (kg.m-2.s-1)
+      do ig=1,ngrid 
+         pdqs_sedco2(ig)=sum_subpdqs_sedco2(ig)/real(imicroco2)
+      enddo
+c Temperature tendency pdtcloud
+      DO l=1,nlay
+        DO ig=1,ngrid
+             pdtcloudco2(ig,l) =
+     &         sum_subpdt(ig,l)/real(imicroco2)-pdt(ig,l)
+        ENDDO
+      ENDDO
+c Tracers tendencies pdqcloud
+      DO l=1,nlay
+        DO ig=1,ngrid        
+             pdqcloudco2(ig,l,igcm_co2_ice) = 
+     &            sum_subpdq(ig,l,igcm_co2_ice)/real(imicroco2)
+     &            - pdq(ig,l,igcm_co2_ice)
+             pdqcloudco2(ig,l,igcm_co2) = 
+     &            sum_subpdq(ig,l,igcm_co2)/real(imicroco2)
+     &            - pdq(ig,l,igcm_co2)
+c D.BARDET :
+             if (co2useh2o) then
+             pdqcloudco2(ig,l,igcm_h2o_ice) = 
+     &            sum_subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
+     &            - pdq(ig,l,igcm_h2o_ice)
+
+             pdqcloudco2(ig,l,igcm_ccn_mass) = 
+     &            sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_mass)
+
+             pdqcloudco2(ig,l,igcm_ccn_number) = 
+     &            sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_number)
+             endif
+        
+             pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
+     &            sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccnco2_mass)
+        
+             pdqcloudco2(ig,l,igcm_ccnco2_number) = 
+     &            sum_subpdq(ig,l,igcm_ccnco2_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccnco2_number)
+
+             pdqcloudco2(ig,l,igcm_dust_mass) = 
+     &            sum_subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_dust_mass)
+
+             pdqcloudco2(ig,l,igcm_dust_number) = 
+     &            sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_dust_number)
+        ENDDO
+      ENDDO
+c Due to stepped entry, other processes tendencies can add up to negative values
+c Therefore, enforce positive values and conserve mass
+      DO l=1,nlay
+        DO ig=1,ngrid
+             IF ((pqeff(ig,l,igcm_ccnco2_number) + 
+     &            ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 
+     &            pdqcloudco2(ig,l,igcm_ccnco2_number))
+     &            .lt. 1.)
+     &            .or. (pqeff(ig,l,igcm_ccnco2_mass) + 
+     &            ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 
+     &            pdqcloudco2(ig,l,igcm_ccnco2_mass))
+     &            .lt. 1.e-20)) THEN
+                pdqcloudco2(ig,l,igcm_ccnco2_number) =
+     &               - pqeff(ig,l,igcm_ccnco2_number)/ptimestep 
+     &               - pdq(ig,l,igcm_ccnco2_number)+1.
+                
+                pdqcloudco2(ig,l,igcm_dust_number) =  
+     &               -pdqcloudco2(ig,l,igcm_ccnco2_number)
+
+                pdqcloudco2(ig,l,igcm_ccnco2_mass) =
+     &               - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
+     &               - pdq(ig,l,igcm_ccnco2_mass)+1.e-20
+
+                pdqcloudco2(ig,l,igcm_dust_mass) = 
+     &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
+             ENDIF
+        ENDDO
+      ENDDO
+      DO l=1,nlay
+        DO ig=1,ngrid
+             IF ( (pqeff(ig,l,igcm_dust_number) + 
+     &            ptimestep* (pdq(ig,l,igcm_dust_number) + 
+     &            pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.)
+     &            .or. (pqeff(ig,l,igcm_dust_mass)+ 
+     &            ptimestep* (pdq(ig,l,igcm_dust_mass) + 
+     &            pdqcloudco2(ig,l,igcm_dust_mass))
+     &            .le. 1.e-20)) then                  
+                pdqcloudco2(ig,l,igcm_dust_number) =
+     &               - pqeff(ig,l,igcm_dust_number)/ptimestep 
+     &               - pdq(ig,l,igcm_dust_number)+1.
+
+                pdqcloudco2(ig,l,igcm_ccnco2_number) =  
+     &               -pdqcloudco2(ig,l,igcm_dust_number)
+
+                pdqcloudco2(ig,l,igcm_dust_mass) =
+     &               - pqeff(ig,l,igcm_dust_mass)/ptimestep
+     &               - pdq(ig,l,igcm_dust_mass) +1.e-20
+
+                pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
+     &               -pdqcloudco2(ig,l,igcm_dust_mass)
+             ENDIF
+        ENDDO
+      ENDDO
+! pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq      
+      DO l=1,nlay
+        DO ig=1,ngrid
+             IF (pqeff(ig,l,igcm_co2_ice) + ptimestep*
+     &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 
+     &       .lt. 1.e-15) THEN
+           pdqcloudco2(ig,l,igcm_co2_ice) = 
+     &     - pqeff(ig,l,igcm_co2_ice)/ptimestep-pdq(ig,l,igcm_co2_ice)
+           pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
+          ENDIF    
+          IF (pqeff(ig,l,igcm_co2) + ptimestep*
+     &       (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 
+     &       .lt. 0.1) THEN
+           pdqcloudco2(ig,l,igcm_co2) = 
+     &     - pqeff(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2)
+           pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2)
+          ENDIF
+        ENDDO
+      ENDDO
+
+c Update clouds parameters values in the cloud fraction (for output)
+      DO l=1, nlay
+        DO ig=1,ngrid
+
+              Niceco2=pqeff(ig,l,igcm_co2_ice) +                   
+     &             (pdq(ig,l,igcm_co2_ice) + 
+     &             pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
+              Nco2=pqeff(ig,l,igcm_co2) +                   
+     &             (pdq(ig,l,igcm_co2) + 
+     &             pdqcloudco2(ig,l,igcm_co2))*ptimestep
+              Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) +                 
+     &             (pdq(ig,l,igcm_ccnco2_number) +               
+     &             pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)
+     &             ,1.e-30)
+              Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) +                 
+     &             (pdq(ig,l,igcm_ccnco2_mass) +               
+     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)
+     &             ,1.e-30)
+              
+              myT=pteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
+              rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 
+     &             myT-2.87e-6* myT* myT)
+              rho_ice_co2=rho_ice_co2T(ig,l)
+c             rho_ice_co2 is shared by tracer_mod and used in updaterice
+c     Compute particle size
+              call updaterice_microco2(Niceco2,
+     &             Qccnco2,Nccnco2,
+     &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
+              
+              if ( (Niceco2 .le. 1.e-25 .or. 
+     &             Nccnco2*tauscaling(ig) .le. 1.) )THEN
+                 riceco2(ig,l)=0.
+                 Qext1bins2(ig,l)=0.
+              else
+c     Compute opacities
+                No=Nccnco2*tauscaling(ig)
+                Rn=-dlog(riceco2(ig,l))
+                n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
+                Qext1bins2(ig,l)=0.
+                do i = 1, nbinco2_cld 
+                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
+                 n_derf = erf((rb_cldco2(i+1)+Rn) *dev2)
+                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
+                 Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
+                enddo
+              endif
+     
+c D.BARDET : update rice water only if co2 use h2o ice as CCN
+          if (co2useh2o) then
+             call updaterice_micro(
+     &       pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &       (pdq(ig,l,igcm_h2o_ice) +                     ! ice mass
+     &       pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &       pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
+     &       (pdq(ig,l,igcm_ccn_mass) +                    ! ccn mass
+     &       pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
+     &       pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
+     &       (pdq(ig,l,igcm_ccn_number) +                  ! ccn number
+     &       pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
+     &       tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+          endif
+
+          call updaterdust(
+     &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
+     &   (pdq(ig,l,igcm_dust_mass) +                     ! dust mass
+     &    pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
+     &    pqeff(ig,l,igcm_dust_number) +                 ! dust number
+     &   (pdq(ig,l,igcm_dust_number) +                   ! dust number
+     &    pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
+     &    rdust(ig,l))
+
+        ENDDO
+      ENDDO 
+     
+c ======================================================================     
+c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
+c Then that should not affect the ice particle radius
+c ======================================================================         
+      do ig=1,ngrid
+        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
+             if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
+     &       riceco2(ig,2)=riceco2(ig,3) 
+             riceco2(ig,1)=riceco2(ig,2)
+        endif
+      end do
+       
+      DO l=1,nlay
+         DO ig=1,ngrid
+           rsedcloud(ig,l)=max(rice(ig,l)*
+     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
+     &                    rdust(ig,l))
+!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
+         ENDDO
+      ENDDO
+       
+      DO l=1,nlay
+         DO ig=1,ngrid
+           rsedcloudco2(ig,l)=max(riceco2(ig,l)*
+     &           (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
+     &                    rdust(ig,l))
+c          rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
+         ENDDO
+      ENDDO
+       
+      call co2sat(ngrid*nlay,pteff+(pdt+pdtcloudco2)*ptimestep
+     &      ,pplay,zqsatco2)
+      do l=1,nlay
+        do ig=1,ngrid 
+             satuco2(ig,l) = (pqeff(ig,l,igcm_co2)  + 
+     &            (pdq(ig,l,igcm_co2) +
+     &            pdqcloudco2(ig,l,igcm_co2))*ptimestep)*
+     &            (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
+        enddo
+      enddo
+! Everything modified by CO2 microphysics must be wrt co2cloudfrac
+      IF (CLFvaryingCO2) THEN
+        DO l=1,nlay
+          DO ig=1,ngrid
+
+            pdqcloudco2(ig,l,igcm_ccnco2_mass)=
+     &        pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_ccnco2_number)=
+     &        pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_dust_mass)=
+     &        pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_dust_number)=
+     &        pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
+c D.BARDET
+            if (co2useh2o) then
+            pdqcloudco2(ig,l,igcm_h2o_ice)=
+     &        pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_ccn_mass)=
+     &        pdqcloudco2(ig,l,igcm_ccn_mass)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_ccn_number)=
+     &        pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)            
+            endif
+
+            pdqcloudco2(ig,l,igcm_co2_ice)=
+     &        pdqcloudco2(ig,l,igcm_co2_ice)*co2cloudfrac(ig,l)
+
+            pdqcloudco2(ig,l,igcm_co2)=
+     &        pdqcloudco2(ig,l,igcm_co2)*co2cloudfrac(ig,l)
+
+            pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*co2cloudfrac(ig,l)
+
+            Qext1bins2(ig,l)=Qext1bins2(ig,l)*co2cloudfrac(ig,l)
+          ENDDO
+        ENDDO   
+      ENDIF
+! opacity in mesh ig is the sum over l of Qext1bins2: Is this true ?
+      tau1mic(:)=0.
+      do l=1,nlay
+        do ig=1,ngrid 
+          tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
+        enddo
+      enddo
+!Outputs:
+      call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
+     &        SatIndex)
+      call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
+     &        satuco2)
+      call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
+     &        ,3,riceco2)         
+      call WRITEdiagfi(ngrid,"co2cloudfrac","co2 cloud fraction"
+     &        ," ",3,co2cloudfrac)
+      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
+     &        ,"m",3,rsedcloudco2)
+      call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
+     &        ," ",3,Qext1bins2)
+      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
+     &        ," ",2,tau1mic)
+      call WRITEDIAGFI(ngrid,"mem_Nccn_co2","CCN number used by CO2"
+     &        ,"kg/kg ",3,mem_Nccn_co2) 
+      call WRITEDIAGFI(ngrid,"mem_Mccn_co2","CCN mass used by CO2"
+     &        ,"kg/kg ",3,mem_Mccn_co2) 
+      call WRITEDIAGFI(ngrid,"mem_Mh2o_co2","H2O mass in CO2 crystal"
+     &        ,"kg/kg ",3,mem_Mh2o_co2)         
+c D.BARDET: sensibility test 
+c      call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussiere"
+c     &        ,"part/kg",3,No_dust)
+c      call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussiere"
+c     &        ,"kg/kg ",3,Mo_dust) 
+      END SUBROUTINE co2cloud
+
+c ===================================================================
+c Subroutines used to write variables of memory in start files       
+c ===================================================================
+
+      SUBROUTINE ini_co2cloud(ngrid,nlayer)
+  
+      IMPLICIT NONE
+
+      INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
+      INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
+
+         allocate(mem_Nccn_co2(ngrid,nlayer))
+         allocate(mem_Mccn_co2(ngrid,nlayer))
+         allocate(mem_Mh2o_co2(ngrid,nlayer))
+
+      END SUBROUTINE ini_co2cloud
+c ----------------------------------
+      SUBROUTINE end_co2cloud
+
+      IMPLICIT NONE
+
+         if (allocated(mem_Nccn_co2)) deallocate(mem_Nccn_co2)
+         if (allocated(mem_Mccn_co2)) deallocate(mem_Mccn_co2)
+         if (allocated(mem_Mh2o_co2)) deallocate(mem_Mh2o_co2)
+
+      END SUBROUTINE end_co2cloud
+
+      END MODULE co2cloud_mod
Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F90	(revision 2348)
+++ 	(revision )
@@ -1,982 +1,0 @@
-!======================================================================================================================!
-! Module: CO2 clouds formation ========================================================================================!
-!----------------------------------------------------------------------------------------------------------------------!
-! Authors: Joachim Audouard, Constantino Listowski, Anni Määttänen
-! Date: 09/2016
-!----------------------------------------------------------------------------------------------------------------------!
-! Contains subroutines:
-!      - co2cloud: of co2 cloud microphysics
-!
-!      - ini_co2cloud: (only used in phys_state_var_init_mod.F90)
-!
-!      - end_co2cloud: (only used in phys_state_var_init_mod.F90)
-!======================================================================================================================!
-module co2cloud_mod
-
-implicit none
-
-double precision, allocatable, save :: &
-   mem_Mccn_co2(:,:), &! Memory of CCN mass of H2O and dust used by CO2
-   mem_Mh2o_co2(:,:), &! Memory of H2O mass integred into CO2 crystal
-   mem_Nccn_co2(:,:)   ! Memory of CCN number of H2O and dust used by CO2
-
-contains
-!======================================================================================================================!
-! SUBROUTINE: co2cloud ================================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Main subroutine of co2 cloud microphysics
-!----------------------------------------------------------------------------------------------------------------------!
-! Comments:
-!----------
-!     - Adaptation of the water ice clouds scheme (with specific microphysics) of Montmessin, Navarro et al.
-!
-!     - Microphysics subroutine is improvedCO2clouds.F
-!
-!     - There is a time loop specific to cloud formation due to timescales smaller than the GCM integration timestep
-!
-!     - The microphysics time step is a fraction of the physical one
-!
-!     - The CO2 clouds tracers (co2_ice, ccn mass and concentration) are sedimented at each microtimestep. pdqs_sedco2
-!       keeps track of the CO2 flux at the surface
-!
-!     - The subgrid Temperature distribution is modulated (0 or 1) by Spiga et al. (GRL 2012)
-!
-!     - Saturation Index to account for GW propagation or dissipation upwards
-!
-!     - 4D and column opacities are computed using Qext values at 1µm
-!----------------------------------------------------------------------------------------------------------------------!
-! Papers:
-!--------
-!   "Near-pure vapor condensation in the Martian atmosphere: CO2 ice crystal growth", Listowski et al. (2013), JGRE
-!   "Modeling the microphysics of CO2 ice clouds within wave-induced cold pockets in the martian mesosphere", Listowski
-!     et al. (2014), Icarus
-!   "Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice
-!     clouds", Navarro et al. (2014), JGRE
-!   "Martian GCM with complete CO2 clouds microphysics", Audouard et al. (2017), EPSC abstract
-!----------------------------------------------------------------------------------------------------------------------!
-! Algorithm:
-!-----------
-!   0. Firstcall
-!     0.1. Initialization of microtimestep from imicroco2
-!     0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
-!     0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
-!     0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins)
-!     0.5. Save the radius grid of CO2 particles (rb_cldco2)
-!   1. Initialization
-!   2. Compute mass and thickness layers
-!   3. Define the sub-grid cloud (CLFvaryingCO2)
-!     3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
-!       3.1.a. Saturation index CO2
-!       3.1.b. Compute tcond
-!       3.1.c. Compute cloud fraction in cells
-!     3.2. No sub-grid cloud representation (CLFvarying = False)
-!   4. Microphysics of CO2 cloud formation
-!     4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
-!     4.2. Effective tracers quantities in the cloud
-!     4.3. Gravitational sedimentation
-!       4.3.a. Compute cloud density
-!       4.3.b. Save actualized tracer values to compute sedimentation tendancies
-!       4.3.c. Sedimentation of co2 ice
-!       4.3.d. Sedimentation for other tracers
-!       4.3.e. Compute tendencies due to the sedimation process
-!     4.4. Main call to the cloud scheme
-!     4.5. Updating tendencies after cloud scheme
-!   5. Compute final tendencies after time loop
-!   6. Update clouds physical values in the cloud (for output)
-!     6.1. Update density of co2 ice, riceco2 and opacity
-!     6.2. Update rice and rdust
-!   7. Correction if a lot of subliming CO2 fills the 1st layer
-!   8. Compute water cloud sedimentation radius
-!   9. CO2 saturated quantities
-!     9.1 Compute CO2 saturation in layers
-!     9.2 Compute CO2 saturated quantities in layers
-!   10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
-!   11. Compute opacity at 1 micron
-!   12. Write outputs in diagfi.nc
-!======================================================================================================================!
-  subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
-                      nq, tau, tauscaling, rdust, rice, riceco2, nuice, rsedcloud, rhocloud, pzlev, pdqs_sedco2, pdu, &
-                      pu, pcondicea, co2ice)
-
-  use ioipsl_getincom, only: getin
-  use dimradmars_mod,  only: naerkind
-  use comcstfi_h,      only: pi, g, cpp
-  use updaterad,       only: updaterice_microco2, updaterice_micro, updaterdust
-  use conc_mod,        only: mmean, rnew
-  use tracer_mod,      only: igcm_co2, igcm_co2_ice, igcm_dust_mass, igcm_dust_number, igcm_h2o_ice, &
-                             igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, rho_dust, &
-                             nuiceco2_sed, nuiceco2_ref, rho_ice_co2, r3n_q, rho_ice, nuice_sed
-
-  use newsedim_mod,    only: newsedim
-
-  use datafile_mod,    only: datadir
-
-  use improvedCO2clouds_mod, only: improvedCO2clouds, density_co2_ice
-
-  implicit none
-
-  include "callkeys.h"
-  include "microphys.h"
-!----------------------------------------------------------------------------------------------------------------------!
-! VARIABLES DECLARATION
-!----------------------------------------------------------------------------------------------------------------------!
-! Input arguments:
-!-----------------
-  integer, intent(in) ::&
-     ngrid, &! Number of grid points
-     nlay,  &! Number of layers
-     nq      ! Number of tracers
-
-  real, intent(in) :: &
-     ptimestep,           &! Physical time step (s)
-     pplev(ngrid,nlay+1), &! Inter-layer pressures (Pa)
-     pplay(ngrid,nlay),   &! Mid-layer pressures (Pa)
-     pdpsrf(ngrid),       &! Tendency on surface pressure
-     pzlay(ngrid,nlay),   &! Altitude at the middle of the layers
-     pt(ngrid,nlay),      &! Temperature at the middle of the layers (K)
-     pdt(ngrid,nlay),     &! Tendency on temperature from other parametrizations
-     pq(ngrid,nlay,nq),   &! Tracers (kg/kg)
-     pdq(ngrid,nlay,nq),  &! Tendencies before condensation (kg/kg.s-1)
-     tau(ngrid,naerkind), &! Column dust optical depth at each point
-     tauscaling(ngrid),   &! Convertion factor for dust amount
-     pu(ngrid,nlay),      &! Zonal Wind: zu = pu + (pdu * ptimestep)
-     pdu(ngrid,nlay),     &! Tendency of zonal wind before condensation
-     pzlev(ngrid,nlay+1), &! Altitude at the boundaries of the layers
-     nuice(ngrid,nlay),   &! Estimated effective variance of the size distribution
-     co2ice(ngrid)         ! Amount of co2 ice at the surface
-!----------------------------------------------------------------------------------------------------------------------!
-! Output arguments:
-!------------------
-  real, intent(out) :: &
-     rdust(ngrid,nlay),         & ! Dust geometric mean radius (m)
-     rice(ngrid,nlay),          & ! Water Ice mass mean radius (m)
-     rsedcloud(ngrid,nlay),     & ! Water Cloud sedimentation radius
-     rhocloud(ngrid,nlay),      & ! Water Cloud density (kg.m-3)
-     pdqs_sedco2(ngrid),        & ! CO2 flux at the surface
-     pdqcloudco2(ngrid,nlay,nq),& ! Tendency due to CO2 condensation (kg/kg.s-1)
-     pcondicea(ngrid,nlay),     & ! Rate of condensation/sublimation of co2 ice in layers
-     pdtcloudco2(ngrid,nlay),   & ! Tendency on temperature due to latent heat
-     riceco2(ngrid,nlay)          ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
-!----------------------------------------------------------------------------------------------------------------------!
-! Local:
-!-------
-!-----1) Parameters:
-!-------------------
-  integer, parameter :: &
-     uQext = 555,         &! file_qext unit ID
-     var_dim_qext = 10000  ! Exact dimension of radv and qextv1mic from file_qext
-
-  real, parameter :: &
-     mincloud = 0.1,  &! Minimum cloud fraction
-     beta = 0.85,     &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95):
-                       !   beta = 1    for spheres
-                       !   beta = 0.85 for irregular particles
-                       !   beta = 0.5  for disk shaped particles
-     threshold = 1e-30 ! limit value
-
-  double precision, parameter :: &
-     rmin_cld = 1.e-9,  &! Minimum radius (m)
-     rmax_cld = 5.e-6,  &! Maximum radius (m)
-     rbmin_cld = 1.e-10,&! Minimum boundary radius (m)
-     rbmax_cld = 2.e-4, &! Maximum boundary radius (m)
-     Fo = 7.5e-7,       &! for sat index  (J.m-3)
-     lambdaH = 150.e3    ! for sat index (km)
-
-  character(len=23), parameter :: &
-     file_qext = 'optprop_co2ice_1mic.dat' ! File extinction coefficients of CO2 particles
-!----------------------------------------------------------------------------------------------------------------------!
-!-----2) Saved:
-!--------------
-  integer, save :: &
-     imicroco2 ! Time subsampling for coupled water microphysics sedimentation microtimestep timeloop for microphysics:
-               !   if imicroco2 = 1, subpdt is the same as pdt
-  real, save :: &
-     sigma_iceco2, &! Variance of the ice and CCN distributions
-     microtimestep  ! Integration timestep for coupled water microphysics & sedimentation
-
-  double precision, save :: &
-     dev2,                   &! 1. / ( sqrt(2.) * sigma_iceco2 )
-     Qext1bins(nbinco2_cld), &! Extinction coefficients for rb_cldco2 radius of CO2 ice particles
-     rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
-
-  logical, save :: &
-     firstcall = .true. ! Used to compute saved variables
-!----------------------------------------------------------------------------------------------------------------------!
-!-----3) Variables:
-!------------------
-  integer :: &
-     iq,       &! loop on tracers
-     ig,       &! loop on grid points
-     l,        &! loop on layers
-     i,        &! loop on nbinco2_cld
-     nelem,    &! number of point between lebon1 and lebon2 => interpolation
-     lebon1,   &! bound limit for the interpolation
-     lebon2,   &! bound limit for the interpolation
-     microstep  ! Time subsampling step variable
-
-  real :: &
-! ---Tendency given by clouds inside the micro loop
-     subpdqcloudco2(ngrid,nlay,nq), &! On tracers, cf. pdqcloud
-     subpdtcloudco2(ngrid,nlay),    &! On temperature, cf. pdtcloud
-!  ---Global tendency (clouds+physics)
-     sum_subpdq(ngrid,nlay,nq),     &! On tracers, cf. pdqcloud
-     sum_subpdt(ngrid,nlay),        &! On temperature, cf. pdtcloud
-!  ---Sedimentation
-     rsedcloudco2(ngrid,nlay),      &! Cloud sedimentation radius
-     ztsed(ngrid,nlay),             &! Temperature with real-time value in microtimeloop
-     zqsed(ngrid,nlay,nq),          &! Tracers with real-time value in µloop
-     zqsed0(ngrid,nlay,nq),         &! For sedimentation tendancy
-     subpdqsed(ngrid,nlay,nq),      &! Tendancy due to sedimentation
-     sum_subpdqs_sedco2(ngrid),     &! CO2 flux at the surface
-!  ---For sub grid T distribution
-     zt(ngrid,nlay),                &! Local value of temperature
-     zq_co2vap(ngrid, nlay),        &! Local value of CO2 vap
-     rhocloudco2t(ngrid, nlay),     &! Cloud density (kg.m-3)
-!  ---For Saturation Index computation
-     zdelt,                         &! Delta T for the temperature distribution
-     co2cloudfrac(ngrid,nlay),      &! Cloud fraction used only with CLFvarying is true
-!  ---Misc
-     rhocloudco2(ngrid, nlay),      &! Cloud density (kg.m-3)
-     Nccnco2,                       &! buffer: number of ccn used for co2 condensation
-     Qccnco2,                       &! buffer: mass of ccn used for co2 condensation
-     Niceco2,                       &! buffer: mmr co2 ice
-     epaisseur(ngrid,nlay),         &! Layer thickness (m)
-     masse(ngrid,nlay),             &! Layer  mass (kg.m-2)
-     pteff(ngrid, nlay),            &! Effective temperature in the cloud
-     pqeff(ngrid, nlay, nq),        &! Effective tracers quantities in the cloud
-     wq(ngrid,nlay+1),              &! Displaced tracer mass (kg.m-2) during microtimestep
-     satuco2(ngrid,nlay),           &! CO2 satu ratio for output diagfi
-     zqsatco2(ngrid,nlay)            ! Saturation co2
-
-  double precision :: &
-! ---Extinction coefficients at 1 micron of CO2 particles
-     vrat_cld,                      &! Volume ratio
-     n_derf,                        &! derf( (rb_cldco2(1)-dlog(riceco2(ig,l))) *dev2)
-     Qtemp,                         &! mean value in the interval during the interpolation
-     ltemp1(var_dim_qext),          &! abs(radv(:)-rb_cldco2(i))
-     ltemp2(var_dim_qext),          &! abs(radv(:)-rb_cldco2(i+1))
-     n_aer(nbinco2_cld),            &! -0.5 * Nccnco2*tauscaling(ig) * n_derf
-     tau1mic(ngrid),                &! CO2 ice column opacity at 1µm
-     Qext1bins2(ngrid,nlay),        &! CO2 ice opacities
-     radv(var_dim_qext),            &! radius of CO2 ice at 1 µm (read from file_qext)
-     Qextv1mic(var_dim_qext),       &! extinction coefficient of CO2 ice at 1 µm (read from file_qext)
-! ---For Saturation Index computation
-     rho,                           &! background density
-     zu,                            &! absolute value of zonal wind field
-     NN,                            &! N^2 static stability
-     gradT,                         &! thermal gradient
-     SatIndex(ngrid,nlay),          &! Saturation index S in Spiga 2012 paper, assuming like in the paper GW phase speed
-                                     !  (stationary waves): c = 0 m.s-1, lambdaH = 150 km, Fo = 7.5e-7 J.m-3
-     SatIndexmap(ngrid),            &! maxval(SatIndex(ig,12:26))
-! ---Misc
-     myT,                           &! Temperature scalar for co2 density computation
-     tcond(ngrid,nlay),             &! CO2 condensation temperature
-     rho_ice_co2T(ngrid,nlay)        ! T-dependant CO2 ice density
-
-  logical :: &
-     file_qext_ok ! Check if file_qext exists
-!======================================================================================================================!
-! BEGIN ===============================================================================================================!
-!======================================================================================================================!
-! 0. Firstcall
-!----------------------------------------------------------------------------------------------------------------------!
-  if (firstcall) then
-    firstcall=.false.
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.1. Initialization of microtimestep from imicroco2
-!----------------------------------------------------------------------------------------------------------------------!
-#ifdef MESOSCALE
-    imicroco2 = 2
-#else
-    imicroco2 = 30
-#endif
-    call getin("imicroco2", imicroco2)
-    microtimestep = ptimestep/real(imicroco2)
-    sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
-    dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
-!        > the grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
-!   - rad_cldco2 is the primary radius grid used for microphysics computation.
-!   - The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
-!   - vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
-!   - The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
-!----------------------------------------------------------------------------------------------------------------------!
-    ! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
-    vrat_cld = exp(log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) * 3.)
-
-    ! rad_cldco2 is the primary radius grid used for microphysics computation.
-    rad_cldco2(1) = rmin_cld
-    do i = 1, nbinco2_cld-1
-      rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
-    end do
-
-    ! rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
-    rb_cldco2(1) = rbmin_cld
-    do i = 1, nbinco2_cld
-      rb_cldco2(i+1) = ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cldco2(i)
-    end do
-    rb_cldco2(nbinco2_cld+1) = rbmax_cld
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
-!----------------------------------------------------------------------------------------------------------------------!
-    ! get information about file_qext
-    inquire(file=trim(datadir)//'/'//file_qext, exist=file_qext_ok)
-
-    ! if file_qext is missing then stop
-    if (.not. file_qext_ok) then
-      write(*,*)'file'//file_qext//'should be in ', trim(datadir)
-      call abort_physic('co2cloud', 'file missing', 1)
-    end if
-
-    ! read file_qext
-    open(unit=uQext,file=trim(datadir)//'/'//file_qext, form='formatted')
-
-    ! skip 1 line
-    read(uQext,*)
-
-    ! extract radv
-    do i = 1, var_dim_qext
-      read(uQext,'(E12.5)')radv(i)
-    end do
-
-    ! skip 1 line
-    read(uQext,*)
-
-    ! Qextv1mic
-    do i = 1 , var_dim_qext
-      read(uQext,'(E12.5)')Qextv1mic(i)
-    end do
-
-    ! close file_qext
-    close(uQext)
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins), using
-!       file_qext values (radv, Qextv1mic)
-!----------------------------------------------------------------------------------------------------------------------!
-    do i = 1, nbinco2_cld
-      ltemp1 = abs(radv(:)-rb_cldco2(i))
-      ltemp2 = abs(radv(:)-rb_cldco2(i+1))
-      lebon1 = minloc(ltemp1,DIM=1)
-      lebon2 = min(minloc(ltemp2,DIM=1), var_dim_qext)
-      nelem = lebon2 - lebon1 + 1.
-
-      ! mean value in the interval
-      Qtemp = 0d0
-      do l = 0, nelem
-        Qtemp = Qtemp + Qextv1mic(min(lebon1+l, var_dim_qext))
-      end do
-
-      Qext1bins(i) = Qtemp / nelem
-    end do
-
-    Qext1bins(:) = Qext1bins(:) * pi * (rad_cldco2(:)**2)
-
-    ! print result of the interpolation
-    write(*,*)'--------------------------------------------'
-    write(*,*)'Microphysics co2: size bin-Qext information:'
-    write(*,*)'   i, rad_cldco2(i), Qext1bins(i)'
-    do i = 1, nbinco2_cld
-      write(*,'(i3,3x,3(e13.6,4x))')i, rad_cldco2(i), Qext1bins(i)
-    end do
-    write(*,*)'--------------------------------------------'
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.5. Save the radius grid of CO2 particles (rb_cldco2)
-!----------------------------------------------------------------------------------------------------------------------!
-    do i = 1, nbinco2_cld+1
-      rb_cldco2(i) = log(rb_cldco2(i))
-    end do
-  end if ! of IF (firstcall)
-!----------------------------------------------------------------------------------------------------------------------!
-! 1. Initialization
-!----------------------------------------------------------------------------------------------------------------------!
-  sum_subpdq(1:ngrid,1:nlay,1:nq) = 0.
-  sum_subpdt(1:ngrid,1:nlay) = 0.
-
-  subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
-  subpdtcloudco2(1:ngrid,1:nlay) = 0.
-
-  pdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
-  pdtcloudco2(1:ngrid,1:nlay) = 0.
-
-  ! default value if no ice
-  rhocloudco2(1:ngrid,1:nlay) = rho_dust
-  rhocloudco2t(1:ngrid,1:nlay) = rho_dust
-  epaisseur(1:ngrid,1:nlay) = 0.
-  masse(1:ngrid,1:nlay) = 0.
-
-  zqsed0(1:ngrid,1:nlay,1:nq) = 0.
-  sum_subpdqs_sedco2(1:ngrid) = 0.
-  subpdqsed(1:ngrid,1:nlay,1:nq) = 0.
-!----------------------------------------------------------------------------------------------------------------------!
-! 2. Compute mass and thickness layers
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig = 1, ngrid
-      masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
-      epaisseur(ig,l) = pzlev(ig,l+1) - pzlev(ig,l)
-    end do
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 3. Define the sub-grid cloud (CLFvaryingCO2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
-!----------------------------------------------------------------------------------------------------------------------!
-  if (CLFvaryingCO2) then
-    ! effective temperature
-    pteff(:,:) = pt(:,:)
-
-    ! min co2cloudfrac when there is ice
-    co2cloudfrac(:,:) = mincloud
-
-    ! temperature
-    do l=1,nlay
-      do ig=1,ngrid
-        zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
-      end do
-    end do
-
-    ! Quantities of traceurs
-    if (igcm_co2 /= 0) then
-      do l = 1, nlay
-        do ig = 1, ngrid
-          zq_co2vap(ig,l) = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep
-        end do
-      end do
-    end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 3.1.a. Saturation index CO2
-!----------------------------------------------------------------------------------------------------------------------!
-    ! if saturation index co2 is true
-    if (satindexco2) then
-      ! layers 12 --> 26 ~ 12->85 km
-      do l = 12, 26
-        do ig = 1, ngrid
-          ! compute N^2 static stability
-          gradT = (zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
-          NN = sqrt(g/zt(iq,l) * (max(gradT,-g/cpp) + g/cpp))
-
-          ! compute absolute value of zonal wind field
-          zu = abs(pu(ig,l) + pdu(ig,l)*ptimestep)
-
-          ! compute background density
-          rho = pplay(ig,l) / (rnew(ig,l)*zt(ig,l))
-
-          ! saturation index: Modulate the DeltaT by GW propagation index:
-          ! --------------------------------------------------------------
-          SatIndex(ig,l) = sqrt(Fo*lambdaH/(2.*pi)*NN / (rho*zu**3) )
-        end do
-      end do
-
-      ! Then compute Satindex map in layers 12 --> 26 ~ 12->85 km
-      do ig = 1, ngrid
-        SatIndexmap(ig) = maxval(SatIndex(ig,12:26))
-      end do
-
-      ! Write outputs in diagfi.nc
-      call WRITEDIAGFI(ngrid, "SatIndex", "SatIndex", " ", 3, SatIndex)
-
-      call WRITEDIAGFI(ngrid, "SatIndexmap", "SatIndexmap", "km", 2, SatIndexmap)
-    !------------------------------------------------------------------------------------------------------------------!
-    ! if saturation index co2 is false, set saturation index to 0.05
-    !------------------------------------------------------------------------------------------------------------------!
-    else
-      do ig = 1, ngrid
-        SatIndexmap(ig)=0.05
-      end do
-    end if ! of if (satindexco2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 3.1.b. Compute tcond
-!----------------------------------------------------------------------------------------------------------------------!
-    call tcondco2(ngrid,nlay,pplay,zq_co2vap,tcond)
-!----------------------------------------------------------------------------------------------------------------------!
-! 3.1.c. Compute cloud fraction in cells
-!----------------------------------------------------------------------------------------------------------------------!
-    do ig = 1, ngrid
-      if (SatIndexmap(ig) <= 0.1) then
-        do l = 1, nlay-1
-
-          ! The entire fraction is saturated
-          if (tcond(ig,l) >= (zt(ig,l)+zdelt) .or. tcond(ig,l) <= 0.) then
-            pteff(ig,l) = zt(ig,l)
-            co2cloudfrac(ig,l) = 1.
-
-          ! No saturation at all
-          else if (tcond(ig,l) <= (zt(ig,l)-zdelt)) then
-            pteff(ig,l) = zt(ig,l) - zdelt
-            co2cloudfrac(ig,l) = mincloud
-
-          ! Mean temperature of the cloud fraction
-          else
-            pteff(ig,l) = (tcond(ig,l)+zt(ig,l)-zdelt) / 2.
-            co2cloudfrac(ig,l) = (tcond(ig,l)-zt(ig,l)+zdelt) / (2.0*zdelt)
-          end if
-
-          pteff(ig,l) = pteff(ig,l) - pdt(ig,l)*ptimestep
-
-          ! check boundary values of co2cloudfrac
-          if (co2cloudfrac(ig,l) <= mincloud) then
-            co2cloudfrac(ig,l) = mincloud
-          else if (co2cloudfrac(ig,l)> 1) then
-            co2cloudfrac(ig,l) = 1.
-          end if
-        end do
-
-      ! SatIndex not favorable for GW: leave pt untouched
-      else
-        pteff(ig,l) = pt(ig,l)
-        co2cloudfrac(ig,l) = mincloud
-      end if ! of if (SatIndexmap <= 0.1)
-    end do ! of ngrid
-    ! TODO: Totalcloud frac of the column missing here
-!----------------------------------------------------------------------------------------------------------------------!
-! 3.2. No sub-grid cloud representation (CLFvarying = False)
-!----------------------------------------------------------------------------------------------------------------------!
-  else
-    do l = 1, nlay
-      do ig = 1, ngrid
-        pteff(ig,l) = pt(ig,l)
-      end do
-    end do
-  end if ! end if (CLFvaryingco2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4. Microphysics of CO2 cloud formation
-!----------------------------------------------------------------------------------------------------------------------!
-  do microstep = 1, imicroco2
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
-!----------------------------------------------------------------------------------------------------------------------!
-    do l = 1, nlay
-      do ig = 1, ngrid
-        ! on temperature
-        sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l)
-
-        ! on tracers
-        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)
-
-        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)
-
-        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass)
-
-        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number)
-
-        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice)
-
-        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)
-
-        if (co2useh2o) then
-          sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + pdq(ig,l,igcm_h2o_ice)
-
-          sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)
-
-          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)
-        end if
-      end do ! ngrid
-    end do ! nlay
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.2. Effective tracers quantities in the cloud
-!----------------------------------------------------------------------------------------------------------------------!
-    pqeff(:,:,:) = pq(:,:,:)
-
-    if (CLFvaryingCO2) then
-      pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:)
-
-      pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:)
-
-      pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:)
-    end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3. Gravitational sedimentation
-!----------------------------------------------------------------------------------------------------------------------!
-    if (sedimentation) then
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3.a. Compute cloud density
-!----------------------------------------------------------------------------------------------------------------------!
-      do l = 1, nlay
-        do ig = 1, ngrid
-          ! temperature during the sedimentation process
-          ztsed(ig,l) = pteff(ig,l) + sum_subpdt(ig,l) * microtimestep
-
-          ! quantities tracers during the sedimentation process
-          zqsed(ig,l,:) = pqeff(ig,l,:) + sum_subpdq(ig,l,:) * microtimestep
-
-          ! density of co2 ice
-          call density_co2_ice(ztsed(ig,l), rho_ice_co2T(ig,l))
-
-          ! assure positive value of co2_ice mmr, ccnco2 number, ccnco2 mass
-          Niceco2 = max(zqsed(ig,l,igcm_co2_ice), threshold)
-          Nccnco2 = max(zqsed(ig,l,igcm_ccnco2_number), threshold)
-          Qccnco2 = max(zqsed(ig,l,igcm_ccnco2_mass), threshold)
-
-          ! Get density cloud and co2 ice particle radius
-          call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
-
-          ! conflict with updaterice_microco2 C.M.
-          if (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1) then
-            riceco2(ig,l) = 1.e-9
-          end if
-
-          ! done in updaterice_microco2 C.M.
-          rhocloudco2t(ig,l) = min( max(rhocloudco2t(ig,l),rho_ice_co2T(ig,l)), rho_dust )
-        end do ! ngrid
-      end do ! nlay
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3.b. Save actualized tracer values to compute sedimentation tendancies
-!----------------------------------------------------------------------------------------------------------------------!
-      zqsed0(:,:,igcm_co2_ice) = zqsed(:,:,igcm_co2_ice)
-      zqsed0(:,:,igcm_ccnco2_mass) = zqsed(:,:,igcm_ccnco2_mass)
-      zqsed0(:,:,igcm_ccnco2_number) = zqsed(:,:,igcm_ccnco2_number)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3.c. Sedimentation of co2 ice
-!----------------------------------------------------------------------------------------------------------------------!
-      rsedcloudco2(1:ngrid, 1:nlay) = 5e-6
-      wq(:,:) = 0.
-      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
-                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_co2_ice), wq, beta)
-
-      do ig = 1, ngrid
-        sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep
-      end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3.d. Sedimentation for other tracers
-!----------------------------------------------------------------------------------------------------------------------!
-      wq(:,:) = 0.
-      ! for ccnco2_mass
-      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, &
-                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_mass), wq, beta)
-      !TODO: ajouter le calcule de la tendance a la surface comme co2ice
-
-      wq(:,:) = 0.
-      ! for ccnco2_number
-      call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, &
-                   rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_number), wq, beta)
-
-      !TODO: ajouter le calcule de la tendance a la surface comme co2ice
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3.e. Compute tendencies due to the sedimation process
-!----------------------------------------------------------------------------------------------------------------------!
-      do l = 1, nlay
-        do ig = 1, ngrid
-          subpdqsed(ig,l,igcm_ccnco2_mass) = ( zqsed(ig,l,igcm_ccnco2_mass) - zqsed0(ig,l,igcm_ccnco2_mass)  ) &
-                                             / microtimestep
-
-          subpdqsed(ig,l,igcm_ccnco2_number) = ( zqsed(ig,l,igcm_ccnco2_number) - zqsed0(ig,l,igcm_ccnco2_number) )&
-                                               / microtimestep
-
-          subpdqsed(ig,l,igcm_co2_ice) = ( zqsed(ig,l,igcm_co2_ice) - zqsed0(ig,l,igcm_co2_ice) ) / microtimestep
-        end do
-      end do
-
-      ! update subtimestep tendencies with sedimentation input
-      do l = 1, nlay
-        do ig = 1, ngrid
-          sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqsed(ig,l,igcm_ccnco2_mass)
-
-          sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + subpdqsed(ig,l,igcm_ccnco2_number)
-
-          sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqsed(ig,l,igcm_co2_ice)
-        end do
-      end do
-    end if !(end if sedimentation)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.4. Main call to the cloud scheme
-!----------------------------------------------------------------------------------------------------------------------!
-    call 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, &
-                           rb_cldco2, sigma_iceco2, dev2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.5. Updating tendencies after cloud scheme
-!----------------------------------------------------------------------------------------------------------------------!
-    do l = 1, nlay
-      do ig = 1, ngrid
-        sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
-
-        sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
-
-        sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
-
-        sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
-
-        sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
-                                              subpdqcloudco2(ig,l,igcm_ccnco2_number)
-
-        sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
-
-        sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
-
-        if (co2useh2o) then
-          sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + subpdqcloudco2(ig,l,igcm_h2o_ice)
-
-          sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + subpdqcloudco2(ig,l,igcm_ccn_mass)
-
-          sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
-        end if
-      end do ! ngrid
-    end do ! nlay
-  end do ! of do microstep = 1, imicroco2
-!----------------------------------------------------------------------------------------------------------------------!
-! 5. Compute final tendencies after time loop
-!----------------------------------------------------------------------------------------------------------------------!
-! condensation/sublimation rate in the atmosphere
-  do l = nlay, 1, -1
-    do ig = 1, ngrid
-      pcondicea(ig,l) = sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2)
-    end do
-  end do
-
-  ! CO2 flux at surface (kg.m-2.s-1)
-  do ig = 1, ngrid
-    pdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) / real(imicroco2)
-  end do
-
-  ! temperature tendency
-  do l = 1, nlay
-    do ig = 1, ngrid
-      pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l)
-    end do
-  end do
-
-  ! tracers tendencies
-  do l = 1, nlay
-    do ig = 1, ngrid
-      pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens_mod4micro
-
-      pdqcloudco2(ig,l,igcm_co2_ice) = ( sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) ) - pdq(ig,l,igcm_co2_ice)
-
-      pdqcloudco2(ig,l,igcm_ccnco2_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2) ) - &
-                                           pdq(ig,l,igcm_ccnco2_mass)
-
-      pdqcloudco2(ig,l,igcm_ccnco2_number) = ( sum_subpdq(ig,l,igcm_ccnco2_number) / real(imicroco2) ) - &
-                                             pdq(ig,l,igcm_ccnco2_number)
-
-      pdqcloudco2(ig,l,igcm_dust_mass) = ( sum_subpdq(ig,l,igcm_dust_mass) / real(imicroco2) ) - &
-                                         pdq(ig,l,igcm_dust_mass)
-
-      pdqcloudco2(ig,l,igcm_dust_number) = ( sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2) ) - &
-                                           pdq(ig,l,igcm_dust_number)
-
-      if (co2useh2o) then
-        pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) -&
-                                         pdq(ig,l,igcm_h2o_ice)
-
-        pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - &
-                                          pdq(ig,l,igcm_ccn_mass)
-
-        pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
-                                            pdq(ig,l,igcm_ccn_number)
-      end if
-    end do ! ngrid
-  end do ! nlay
-!----------------------------------------------------------------------------------------------------------------------!
-! 6. Update clouds physical values in the cloud (for output)
-!----------------------------------------------------------------------------------------------------------------------!
-! 6.1. Update density of co2 ice, riceco2 and opacity
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig = 1, ngrid
-      Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep
-
-      Nccnco2 = max( (pqeff(ig,l,igcm_ccnco2_number) + (pdq(ig,l,igcm_ccnco2_number) + &
-                       pdqcloudco2(ig,l, igcm_ccnco2_number)) * ptimestep) &
-                    , threshold)
-
-      Qccnco2 = max( (pqeff(ig,l,igcm_ccnco2_mass) + (pdq(ig,l,igcm_ccnco2_mass) + &
-                       pdqcloudco2(ig,l, igcm_ccnco2_mass)) * ptimestep)&
-                    , threshold)
-
-      myT = pteff(ig,l) + (pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
-
-      ! compute density of co2 ice
-      call density_co2_ice(myT, rho_ice_co2T(ig,l))
-
-      rho_ice_co2 = rho_ice_co2T(ig,l) ! rho_ice_co2 is shared by tracer_mod and used in updaterice
-
-      ! Compute particle size
-      call updaterice_microco2(Niceco2, Qccnco2, Nccnco2, tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
-
-     ! Compute opacities
-      if ( (Niceco2 <= threshold .or. Nccnco2*tauscaling(ig) <= 1.) ) then
-        riceco2(ig,l) = 0.
-        Qext1bins2(ig,l) = 0.
-      else
-        n_derf = derf( (rb_cldco2(1)-dlog(riceco2(ig,l))) *dev2)
-        Qext1bins2(ig,l) = 0.
-
-        do i = 1, nbinco2_cld
-          n_aer(i) = -0.5 * Nccnco2*tauscaling(ig) * n_derf
-
-          n_derf = derf((rb_cldco2(i+1)-dlog(riceco2(ig,l))) *dev2)
-          n_aer(i) = n_aer(i) + (0.5 * Nccnco2*tauscaling(ig) * n_derf)
-
-          Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
-        end do
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 6.2. Update rice and rdust
-!----------------------------------------------------------------------------------------------------------------------!
-      ! update rice water only if co2 use h2o ice as CCN
-      if (co2useh2o) then
-        call updaterice_micro( &
-              pqeff(ig,l,igcm_h2o_ice) + (pdq(ig,l,igcm_h2o_ice) + pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, &
-              pqeff(ig,l,igcm_ccn_mass) + (pdq(ig,l,igcm_ccn_mass) + pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, &
-              pqeff(ig,l,igcm_ccn_number) + (pdq(ig,l,igcm_ccn_number) + pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, &
-              tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-      end if
-
-      ! update rdust
-      call updaterdust( &
-           pqeff(ig,l,igcm_dust_mass) + (pdq(ig,l,igcm_dust_mass) + pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, &
-           pqeff(ig,l,igcm_dust_number) + (pdq(ig,l,igcm_dust_number) + pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, &
-           rdust(ig,l))
-    end do ! ngrid
-  end do ! nlay
-!----------------------------------------------------------------------------------------------------------------------!
-! 7. A correction if a lot of subliming CO2 fills the 1st layer FF (04/2005). Then that should not affect the ice
-!      particle radius
-!----------------------------------------------------------------------------------------------------------------------!
-  do ig = 1, ngrid
-    if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,2))) then
-
-      if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,3)) ) then
-        riceco2(ig,2) = riceco2(ig,3)
-      end if
-
-      riceco2(ig,1) = riceco2(ig,2)
-    end if
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 8. Compute water cloud sedimentation radius
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig= 1, ngrid
-      rsedcloud(ig,l) = max( rice(ig,l)*(1.+nuice_sed)* (1.+nuice_sed)*(1.+nuice_sed), rdust(ig,l) )
-    end do
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 9. CO2 saturated quantities
-!----------------------------------------------------------------------------------------------------------------------!
-! 9.1 Compute CO2 saturation in layers
-!----------------------------------------------------------------------------------------------------------------------!
-  call co2sat(ngrid*nlay, pteff+(pdt+pdtcloudco2)*ptimestep, zqsatco2)
-!----------------------------------------------------------------------------------------------------------------------!
-! 9.2 Compute CO2 saturated quantities in layers
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig = 1, ngrid
-      satuco2(ig,l) = ( pqeff(ig,l,igcm_co2)  + (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))*ptimestep ) *  &
-                      (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) / zqsatco2(ig,l)
-    end do
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
-!----------------------------------------------------------------------------------------------------------------------!
-  if (CLFvaryingCO2) then
-    do l = 1, nlay
-      do ig = 1, ngrid
-        pdqcloudco2(ig,l,igcm_ccnco2_mass) = pdqcloudco2(ig,l,igcm_ccnco2_mass) * co2cloudfrac(ig,l)
-
-        pdqcloudco2(ig,l,igcm_ccnco2_number) = pdqcloudco2(ig,l,igcm_ccnco2_number) * co2cloudfrac(ig,l)
-
-        pdqcloudco2(ig,l,igcm_dust_mass) = pdqcloudco2(ig,l,igcm_dust_mass) * co2cloudfrac(ig,l)
-
-        pdqcloudco2(ig,l,igcm_dust_number) = pdqcloudco2(ig,l,igcm_dust_number) * co2cloudfrac(ig,l)
-
-        pdqcloudco2(ig,l,igcm_co2_ice) = pdqcloudco2(ig,l,igcm_co2_ice) * co2cloudfrac(ig,l)
-
-        pdqcloudco2(ig,l,igcm_co2) = pdqcloudco2(ig,l,igcm_co2) * co2cloudfrac(ig,l)
-
-        pdtcloudco2(ig,l) = pdtcloudco2(ig,l) * co2cloudfrac(ig,l)
-
-        Qext1bins2(ig,l) = Qext1bins2(ig,l) * co2cloudfrac(ig,l)
-
-        if (co2useh2o) then
-          pdqcloudco2(ig,l,igcm_h2o_ice) = pdqcloudco2(ig,l,igcm_h2o_ice) * co2cloudfrac(ig,l)
-
-          pdqcloudco2(ig,l,igcm_ccn_mass) = pdqcloudco2(ig,l,igcm_ccn_mass) * co2cloudfrac(ig,l)
-
-          pdqcloudco2(ig,l,igcm_ccn_number) = pdqcloudco2(ig,l,igcm_ccn_number) * co2cloudfrac(ig,l)
-        end if
-      end do ! ngrid
-    end do ! nlay
-  end if ! if CLFvaryingCO2 is true
-!----------------------------------------------------------------------------------------------------------------------!
-! 11. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
-!----------------------------------------------------------------------------------------------------------------------!
-  tau1mic(:)=0.
-  do l = 1, nlay
-    do ig = 1, ngrid
-      tau1mic(ig) = tau1mic(ig) + Qext1bins2(ig,l)
-    end do
-  end do
-!----------------------------------------------------------------------------------------------------------------------!
-! 12. Write outputs in diagfi.nc
-!----------------------------------------------------------------------------------------------------------------------!
-  call WRITEDIAGFI(ngrid, "satuco2", "vap in satu", " ", 3, satuco2)
-
-  call WRITEDIAGFI(ngrid, "co2cloudfrac", "co2 cloud fraction", " ", 3, co2cloudfrac)
-
-  call writediagfi(ngrid, "riceco2", "ice radius", "m", 3, riceco2)
-
-  call WRITEDIAGFI(ngrid, "Tau3D1mic", " co2 ice opacities", " ", 3, Qext1bins2)
-
-  call WRITEDIAGFI(ngrid, "tau1mic", "co2 ice opacity 1 micron", " ", 2, tau1mic)
-
-  call WRITEDIAGFI(ngrid, "mem_Nccn_co2", "CCN number used by CO2", "kg/kg", 3, mem_Nccn_co2)
-
-  call WRITEDIAGFI(ngrid, "mem_Mccn_co2", "CCN mass used by CO2",  "kg/kg", 3, mem_Mccn_co2)
-
-  call WRITEDIAGFI(ngrid, "mem_Mh2o_co2", "H2O mass in CO2 crystal", "kg/kg", 3, mem_Mh2o_co2)
-!======================================================================================================================!
-! END =================================================================================================================!
-!======================================================================================================================!
-  end subroutine co2cloud
-
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-
-!======================================================================================================================!
-! SUBROUTINE: ini_co2cloud ============================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Allocate arrays used for co2 microphysics
-!======================================================================================================================!
-  subroutine ini_co2cloud(ngrid,nlayer)
-
-  implicit none
-
-  integer, intent(in) :: &
-     ngrid, &! number of atmospheric columns
-     nlayer  ! number of atmospheric layers
-
-  allocate(mem_Nccn_co2(ngrid,nlayer))
-  allocate(mem_Mccn_co2(ngrid,nlayer))
-  allocate(mem_Mh2o_co2(ngrid,nlayer))
-
-  end subroutine ini_co2cloud
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-!======================================================================================================================!
-! SUBROUTINE: end_co2cloud ============================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Deallocate arrays used for co2 microphysics
-!======================================================================================================================!
-  subroutine end_co2cloud
-
-  implicit none
-
-  if (allocated(mem_Nccn_co2)) deallocate(mem_Nccn_co2)
-  if (allocated(mem_Mccn_co2)) deallocate(mem_Mccn_co2)
-  if (allocated(mem_Mh2o_co2)) deallocate(mem_Mh2o_co2)
-
-  end subroutine end_co2cloud
-
-end module co2cloud_mod
Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod4micro.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2condens_mod4micro.F90	(revision 2348)
+++ 	(revision )
@@ -1,685 +1,0 @@
-!======================================================================================================================!
-! Module: CO2 condensation for the CO2 cloud microphysics =============================================================!
-!----------------------------------------------------------------------------------------------------------------------!
-! Authors: Christophe Mathé, Anni Määttänen
-! Date: 16/04/2020
-!----------------------------------------------------------------------------------------------------------------------!
-! Contains subroutines:
-!     - co2condens4micro: condensation/sublimation of CO2 ice on the ground and compute pressure change resulting
-!
-!     - vl1d: Van-Leer scheme
-!======================================================================================================================!
-module co2condens_mod4micro
-
-implicit none
-
-contains
-!======================================================================================================================!
-! SUBROUTINE: co2condens4micro ========================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Condensation/sublimation of CO2 ice on the ground and compute pressure change resulting
-!----------------------------------------------------------------------------------------------------------------------!
-! Comments:
-!----------
-!   Adapted from co2condens_mod.F
-!----------------------------------------------------------------------------------------------------------------------!
-! Paper:
-!-------
-!  Forget et al. (2008), "Non condensable gas enrichment and depletion in the Martian polar regions."
-!----------------------------------------------------------------------------------------------------------------------!
-! Algorithm:
-!-----------
-!   1. Initialization
-!   2. Firstcall
-!   3. Compute CO2 Volume mixing ratio
-!   4. Set zcondicea, zfallice from co2clouds condensation rate and set zt
-!   5. Main co2condens
-!     5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol
-!     5.2. Check if we have condensation/sublimation on the ground
-!     5.3. Compute zfallheat
-!     5.4. Compute direct condensation/sublimation of CO2 ice
-!       5.4.a. If there is not enough CO2 tracer in 1st layer to condense
-!       5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
-!     5.5. Changing CO2 ice amount and pressure
-!     5.6. Surface albedo and emissivity of the surface below the snow (emisref)
-!       5.6.a. Check that amont of CO2 ice is not problematic
-!       5.6.b. Set albedo and emissivity of the surface
-!       5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow)
-!     5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and
-!         warming/cooling of the CO2 currently changing phase).
-!       5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
-!       5.7.b. Mass of each layer at the end of timestep
-!       5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport)
-!         5.7.c.i. Value transfert at the surface interface when condensation/sublimation
-!         5.7.c.ii. Van Leer scheme
-!         5.7.c.iii Compute tendencies on T, U, V, Q
-!   6. CO2 snow / clouds scheme
-!   7. Extra special case for surface temperature tendency pdtsrfc
-!----------------------------------------------------------------------------------------------------------------------!
-  subroutine co2condens4micro(ngrid, nlayer, nq, ptimestep, pcapcal, pplay, pplev, ptsrf, pt, pphi, pdt, pdu, pdv, &
-                              pdtsrf, pu, pv, pq, pdq, piceco2, psolaralb, pemisurf, pdtc, pdtsrfc, pdpsrf, pduc, &
-                              pdvc, pdqc, fluxsurf_sw, zls, zdqssed_co2, pcondicea_co2microp, zdtcloudco2)
-
-  use tracer_mod, only: noms, igcm_co2, igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_number, igcm_dust_mass, &
-                        igcm_ccnco2_number, igcm_ccnco2_mass, igcm_ccn_mass
-
-  use surfdat_h, only: emissiv, phisfi
-
-  use geometry_mod, only: latitude, longitude_deg, latitude_deg
-
-  use planete_h, only: obliquit
-
-  use comcstfi_h, only: cpp, g, r, pi
-
-#ifndef MESOSCALE
-  use vertical_layers_mod, only: ap, bp
-#endif
-
-  implicit none
-
-  include "callkeys.h"
-!----------------------------------------------------------------------------------------------------------------------!
-! VARIABLE DECLARATION
-!----------------------------------------------------------------------------------------------------------------------!
-! Input arguments:
-!-----------------
-  integer, intent(in) :: &
-     nq,    &! number of tracers
-     ngrid, &! number of atmospheric columns
-     nlayer  ! number of vertical layers
-
-  real, intent(in) :: &
-     ptimestep,                       &! physics timestep (s)
-     pcapcal(ngrid),                  &! surface specific heat
-     pplay(ngrid,nlayer),             &! mid-layer pressure (Pa)
-     pplev(ngrid,nlayer+1),           &! inter-layer pressure (Pa)
-     ptsrf(ngrid),                    &! surface temperature (K)
-     pphi(ngrid,nlayer),              &! geopotential (m2.s-2)
-     pt(ngrid,nlayer),                &! atmospheric temperature (K)
-     pu(ngrid,nlayer),                &! zonal wind (m/s)
-     pv(ngrid,nlayer),                &! meridional wind (m/s)
-     pdt(ngrid,nlayer),               &! tendency on temperature from previous physical processes (K/s)
-     pdu(ngrid,nlayer),               &! tendency on zonal wind from previous physical processes(m/s2)
-     pdv(ngrid,nlayer),               &! tendency on meridional wind from previous physical processes (m/s2)
-     pdtsrf(ngrid),                   &! tendency on surface temperature from previous physical processes (K/s)
-     pq(ngrid,nlayer,nq),             &! tracers (../kg_air)
-     pdq(ngrid,nlayer,nq),            &! tendency on tracers from previous physical processes
-     zdqssed_co2(ngrid),              &! CO2 flux at the surface  (kg.m-2.s-1)
-     zdtcloudco2(ngrid,nlayer),       &! tendency on temperature due to latent heat
-     fluxsurf_sw(ngrid,2),            &! added to calculate flux dependent albedo
-     zls,                             &! solar longitude (rad)
-     pcondicea_co2microp(ngrid,nlayer) ! tendency due to CO2 condensation (kg/kg.s-1)
-!----------------------------------------------------------------------------------------------------------------------!
-! In/output arguments:
-!---------------------
-  real, intent(inout) :: &
-     piceco2(ngrid),   &! CO2 ice on the surface (kg.m-2)
-     pemisurf(ngrid),  &! emissivity of the surface
-     psolaralb(ngrid,2) ! albedo of the surface
-!----------------------------------------------------------------------------------------------------------------------!
-! Output arguments:
-!------------------
-  real, intent(out) :: &
-     pdtc(ngrid,nlayer),  &! tendency on temperature dT/dt due to cond/sub (K/s)
-     pdtsrfc(ngrid),      &! tendency on surface temperature (K/s)
-     pdpsrf(ngrid),       &! tendency on surface pressure (Pa/s)
-     pduc(ngrid,nlayer),  &! tendency on zonal wind (m.s-2)
-     pdvc(ngrid,nlayer),  &! tendency on meridional wind (m.s-2)
-     pdqc(ngrid,nlayer,nq) ! tendency on tracers
-!----------------------------------------------------------------------------------------------------------------------!
-! Local:
-!-------
-!----1) Parameters:
-!------------------
-  real, parameter :: &
-     latcond = 595594, &! latent heat of solid CO2 ice (J/kg)
-     cpice = 1000.,    &! specific heat of CO2 ice (J.kg-1.K-1)
-     tcond1mb = 136.27,&! condensation temperature at 1 mbar (K)
-     m_co2 = 44.01E-3, &! CO2 molecular mass (kg/mol)
-     m_noco2 = 33.37E-3 ! non condensible molecular mass (kg/mol)
-
-  logical, parameter :: &
-     improved_ztcond = .true. ! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer) then
-                              !  condensation temperature is computed using partial pressure of CO2
-!----------------------------------------------------------------------------------------------------------------------!
-!----2) Saved:
-!-------------
-  real, save :: &
-     A,    &! coefficient used to compute mean molecular mass
-     B,    &! coefficient used to compute mean molecular mass
-     acond,&! coefficient used to compute ztcondsol
-     bcond,&! coefficient used to compute ztcondsol
-     ccond  ! coefficient used to compute ztcondsol
-
-  logical, save :: &
-     firstcall = .true. ! Used to compute saved variables
-!----------------------------------------------------------------------------------------------------------------------!
-!----3) Variables:
-!-----------------
-  integer :: &
-     l,  &! loop on layers
-     ig, &! loop on ngrid points
-     iq   ! loop on tracer
-
-  real :: &
-     qco2,                   &! effective quantity of CO2, used to compute mean molecular mass
-     mmean,                  &! mean molecular mass
-     zfallheat,              &! aerodynamical friction and energy used to heat the amount of ice
-     zmflux(nlayer+1),       &! mass fluxes through the sigma levels (kg.m-2.s-1)
-     zu(nlayer),             &! effective zonal wind
-     zv(nlayer),             &! effective meridional wind
-     zq(nlayer,nq),          &! effective tracers quantities
-     zq1(nlayer),            &! buffer of zq
-     ztsrf(ngrid),           &! effective temperature at the surface
-     ztm(nlayer+1),          &! temperature fluxes through the sigma levels
-     zum(nlayer+1),          &! zonal wind fluxes through the sigma levels
-     zvm(nlayer+1),          &! meridional wind fluxes through the sigma levels
-     zqm(nlayer+1,nq),       &! quantity of tracers flux through the sigma levels
-     zqm1(nlayer),           &! quantity of tracers after Van-Leer scheme
-     masse(nlayer),          &! mass layer (kg)
-     w(nlayer+1),            &! total mass fluxes through the sigma levels during ptimestep (kg.m-2)
-     availco2,               &! available quantity of co2 (kg)
-     emisref(ngrid),         &! emissivity reference
-     vmr_co2(ngrid,nlayer),  &! CO2 volume mixing ratio
-     zt(nlayer),             &! effective temperature in the atmosphere (K)
-     ztcond(ngrid,nlayer+1), &! CO2 condensation temperature (atm)
-     ztcondsol(ngrid),       &! CO2 condensation temperature (surface)
-     zdiceco2(ngrid),        &! tendency on co2ice surface tracer (kg/m2/s)
-     zcondicea(ngrid,nlayer),&! condensation rate in layer l (kg/m2/s)
-     zcondices(ngrid),       &! condensation rate on the ground (kg/m2/s)
-     zfallice(ngrid)          ! amount of ice falling from layer l (kg/m2/s)
-
-  logical :: &
-     condsub(ngrid) ! True if there is condensation/sublimation (used for co2snow)
-!======================================================================================================================!
-! BEGIN ===============================================================================================================!
-!======================================================================================================================!
-! 1. Initialization
-!----------------------------------------------------------------------------------------------------------------------!
-  availco2 = 0.
-  zfallheat = 0.
-  zt(1:nlayer) = 0.
-  ztcond(1:ngrid, 1:nlayer+1) = 0.
-  ztcondsol(1:ngrid) = 0.
-  zmflux(1:nlayer+1) = 0.
-  zu(1:nlayer) = 0.
-  zv(1:nlayer) = 0.
-  zq(1:nlayer, 1:nq) = 0.
-  zq1(1:nlayer) = 0.
-  ztsrf(1:ngrid) = 0.
-  ztm(1:nlayer+1) = 0.
-  zum(1:nlayer+1) = 0.
-  zvm(1:nlayer+1) = 0.
-  zqm(1:nlayer+1, 1:nq) = 0.
-  masse(1:nlayer) = 0.
-  w(1:nlayer+1) = 0.
-  emisref(1:ngrid) = 0.
-  vmr_co2(1:ngrid, 1:nlayer) = 0.
-  zcondices(1:ngrid) = 0.
-  pdtsrfc(1:ngrid) = 0.
-  pdpsrf(1:ngrid) = 0.
-  zdiceco2(1:ngrid) = 0.
-  condsub(1:ngrid) = .false.
-  zcondicea(1:ngrid, 1:nlayer) = 0.
-  zfallice(1:ngrid) = 0.
-  pduc(1:ngrid, 1:nlayer) = 0.
-  pdvc(1:ngrid, 1:nlayer) = 0.
-  pdqc(1:ngrid, 1:nlayer, 1:nq) = 0.
-  pdtc(1:ngrid,1:nlayer) = 0.
-!----------------------------------------------------------------------------------------------------------------------!
-! 2. Firstcall
-!----------------------------------------------------------------------------------------------------------------------!
-! AS: firstcall OK absolute
-  if (firstcall) then
-    firstcall = .false.
-
-    bcond = 1. / tcond1mb
-    ccond = cpp / (g*latcond)
-    acond = r / latcond
-
-    write(*,*)'CO2condens: improved_ztcond=', improved_ztcond
-    write(*,*)'In co2condens:Tcond(P=1mb)=', tcond1mb, ' Lcond=', latcond
-    write(*,*)'acond,bcond,ccond', acond, bcond, ccond
-
-!   Prepare Special treatment if one of the tracer is CO2 gas. Compute A and B coefficient use to compute mean molecular
-!   mass Mair defined by:
-!     1/Mair = q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2
-!     1/Mair = A*q(igcm_co2) + B
-    A = (1./m_co2 - 1./m_noco2)
-    B = 1./m_noco2
-  end if ! of IF (firstcall)
-!----------------------------------------------------------------------------------------------------------------------!
-! 3. Compute CO2 Volume mixing ratio
-!----------------------------------------------------------------------------------------------------------------------!
-  if (improved_ztcond.and.(igcm_co2/=0)) then
-    do l = 1, nlayer
-      do ig = 1, ngrid
-        qco2 = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep
-!       Mean air molecular mass = 1/(q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2)
-        mmean = 1. / (A*qco2 +B)
-        vmr_co2(ig,l) = (qco2*mmean) / m_co2
-      end do
-    end do
-  else
-    do l = 1, nlayer
-      do ig = 1, ngrid
-        vmr_co2(ig,l) = 0.95
-      end do
-    end do
-  end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 4. Set zcondicea, zfallice from co2clouds condensation rate
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = nlayer, 1, -1
-    do ig = 1, ngrid
-      zcondicea(ig,l) = pcondicea_co2microp(ig,l) *  (pplev(ig,l) - pplev(ig,l+1))/g
-    end do
-  end do
-
-! Only sedimentation falls on the ground !
-  do ig = 1, ngrid
-      zfallice(ig) =  zdqssed_co2(ig)
-      piceco2(ig) = piceco2(ig) + zfallice(ig)*ptimestep
-  end do
-
-  call writediagfi(ngrid, "zcondicea","Condensation rate in layers", " ", 3, zcondicea)
-
-  call writediagfi(ngrid,"zfallice", "Sedimentation rate", " ", 2, zfallice)
-!----------------------------------------------------------------------------------------------------------------------!
-! 5. Main co2condens
-!----------------------------------------------------------------------------------------------------------------------!
-  do ig = 1, ngrid
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol
-!----------------------------------------------------------------------------------------------------------------------!
-    ztcondsol(ig) = 1. / (bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1)))
-    ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.2. Check if we have condensation/sublimation on the ground
-!----------------------------------------------------------------------------------------------------------------------!
-!        ground condensation       ||   falling snow       ||    ground sublimation
-!----------------------------------------------------------------------------------------------------------------------!
-    if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zfallice(ig)/=0.) .OR. ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then
-      condsub(ig) = .true.
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.3. Compute zfallheat
-!----------------------------------------------------------------------------------------------------------------------!
-      if (zfallice(ig)>0) then
-        zfallheat = zfallice(ig) * ( pphi(ig,1) - phisfi(ig) + cpice*(ztcond(ig,1) - ztcondsol(ig)) )/latcond
-      else
-        zfallheat = 0.
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.4. Compute direct condensation/sublimation of CO2 ice
-!----------------------------------------------------------------------------------------------------------------------!
-      zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) / (latcond*ptimestep) - zfallheat
-      pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.4.a. If there is not enough CO2 tracer in 1st layer to condense
-!----------------------------------------------------------------------------------------------------------------------!
-      if (igcm_co2/=0) then
-!       Available CO2 tracer in layer 1 at end of timestep (kg/m2)
-        availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) + (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)*ptimestep) )
-        if ( zcondices(ig) * ptimestep>availco2 ) then
-          zcondices(ig) = availco2/ptimestep
-          zdiceco2(ig) = zcondices(ig) + zfallice(ig)
-          pdtsrfc(ig) = (latcond/pcapcal(ig)) *  (zcondices(ig)+zfallheat)
-        end if
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
-!----------------------------------------------------------------------------------------------------------------------!
-      if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then
-        zcondices(ig) = -piceco2(ig)/ptimestep
-        pdtsrfc(ig) = (latcond/pcapcal(ig)) *  (zcondices(ig)+zfallheat)
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.5. Changing CO2 ice amount and pressure
-!----------------------------------------------------------------------------------------------------------------------!
-      zdiceco2(ig) = zcondices(ig) + zfallice(ig) + sum(zcondicea(ig,:))
-      piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
-      pdpsrf(ig) = -zdiceco2(ig) * g
-
-      if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then
-        print *, 'STOP in condens'
-        print *, 'condensing more than total mass'
-        print *, 'Grid point ', ig
-        print *, 'Longitude(degrees): ', longitude_deg(ig)
-        print *, 'Latitude(degrees): ', latitude_deg(ig)
-        print *, 'Ps = ', pplev(ig,1)
-        print *, 'd Ps = ', pdpsrf(ig)
-        call abort_physic('co2condens4micro', 'condensing more than total mass', 1)
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.6. Surface albedo and emissivity of the surface below the snow (emisref)
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.6.a. Check that amont of CO2 ice is not problematic
-!----------------------------------------------------------------------------------------------------------------------!
-      if(.not.piceco2(ig)>=0.) then
-        if(piceco2(ig)<=-5.e-8) then
-          write(*,*)'WARNING co2condens piceco2(', ig, ')=', piceco2(ig)
-          piceco2(ig) = 0.
-        end if
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.6.b. Set albedo and emissivity of the surface
-!----------------------------------------------------------------------------------------------------------------------!
-      call albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow)
-!----------------------------------------------------------------------------------------------------------------------!
-      if (piceco2(ig)==0) then
-        pemisurf(ig) = emissiv
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!  5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and
-!         warming/cooling of the CO2 currently changing phase).
-!----------------------------------------------------------------------------------------------------------------------!
-      do l= 1, nlayer
-        zt(l) = pt(ig,l) + pdt(ig,l)*ptimestep
-        zu(l) = pu(ig,l) + pdu(ig,l)*ptimestep
-        zv(l) = pv(ig,l) + pdv(ig,l)*ptimestep
-        do iq=1,nq
-          zq(l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-        end do
-      end do
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
-!----------------------------------------------------------------------------------------------------------------------!
-      zmflux(1) = - zcondices(ig) - zfallice(ig)
-      do l = 1, nlayer
-        zmflux(l+1) = zmflux(l) - zcondicea(ig,l) &
-#ifndef MESOSCALE
-                      + (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g)
-! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
-        if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then
-          zmflux(l+1) = 0.
-        end if
-#else
-        if (abs(zmflux(l+1))<1E-13) then
-          zmflux(l+1) = 0.
-        end if
-#endif
-      end do
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.b. Mass of each layer at the end of timestep
-!----------------------------------------------------------------------------------------------------------------------!
-      do l = 1, nlayer
-        masse(l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1))*pdpsrf(ig)*ptimestep)/g
-      end do
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport)
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.c.i. Value transfert at the surface interface when condensation/sublimation
-!----------------------------------------------------------------------------------------------------------------------!
-      ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
-      zum(1) = 0
-      zvm(1) = 0
-
-!     Most tracer do not condense
-      do iq = 1, nq
-        zqm(1,iq) = 0.
-      end do
-
-!     Special case if one of the tracer is CO2 gas
-      if (igcm_co2/=0) then
-        zqm(1,igcm_co2) = 1. ! flux is 100% CO2
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.c.ii. Van Leer scheme
-!----------------------------------------------------------------------------------------------------------------------!
-      do l=1,nlayer
-        w(l)=-zmflux(l)*ptimestep
-      end do
-
-      call vl1d(nlayer,zt,2.,masse,w,ztm)
-      call vl1d(nlayer,zu ,2.,masse,w,zum)
-      call vl1d(nlayer,zv ,2.,masse,w,zvm)
-
-      do iq=1, nq
-        do l=1, nlayer
-          zq1(l) = zq(l,iq)
-        end do
-
-        zqm1(1) = zqm(1,iq)
-        zqm1(2:nlayer) = 0.
-
-        call vl1d(nlayer,zq1,2.,masse,w,zqm1)
-
-        do l = 2, nlayer
-          zq(l,iq) = zq1(l)
-          zqm(l,iq) = zqm1(l)
-        end do
-      end do
-
-!     Surface condensation affects low winds
-      if (zmflux(1)<0) then
-        zum(1) = zu(1) * (w(1)/masse(1))
-        zvm(1) = zv(1) * (w(1)/masse(1))
-        if (w(1)>masse(1)) then ! ensure numerical stability
-          zum(1) = ((zu(1)-zum(2))*masse(1)/w(1)) + zum(2)
-          zvm(1) = ((zv(1)-zvm(2))*masse(1)/w(1)) + zvm(2)
-        end if
-      end if
-
-      ztm(nlayer+1) = zt(nlayer) ! should not be used, but...
-      zum(nlayer+1) = zu(nlayer)  ! should not be used, but...
-      zvm(nlayer+1) = zv(nlayer)  ! should not be used, but...
-
-      do iq = 1, nq
-        zqm(nlayer+1,iq) = zq(nlayer,iq)
-      end do
-!----------------------------------------------------------------------------------------------------------------------!
-!   5.7.c.iii Compute tendencies on T, U, V, Q
-!----------------------------------------------------------------------------------------------------------------------!
-#ifdef MESOSCALE
-! AS: This part must be commented in the mesoscale model
-! AS: ... to avoid instabilities.
-! AS: you have to compile with -DMESOSCALE to do so
-#else
-      do l = 1, nlayer
-        if (.not. co2clouds) then
-!         Tendencies on T
-          pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l))  - zmflux(l+1)*(ztm(l+1) - zt(l)) &
-                         + zcondicea(ig,l)*(ztcond(ig,l)-zt(l)) )
-        else
-          pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l))  - zmflux(l+1)*(ztm(l+1) - zt(l)) )
-        end if
-
-!       Tendencies on U
-        pduc(ig,l) = (1./masse(l)) * ( zmflux(l)*(zum(l) - zu(l)) - zmflux(l+1)*(zum(l+1) - zu(l)) )
-
-!       Tendencies on V
-        pdvc(ig,l) = (1./masse(l)) * ( zmflux(l)*(zvm(l) - zv(l)) - zmflux(l+1)*(zvm(l+1) - zv(l)) )
-      end do
-#endif
-!       Tendencies on Q
-      do iq = 1, nq
-        if (iq==igcm_co2) then
-          do l = 1, nlayer
-            pdqc(ig,l,iq) = (1./masse(l)) * (zmflux(l)*(zqm(l,iq) - zq(l,iq))- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
-                            + zcondicea(ig,l)*(zq(l,iq) - 1.))
-          end do
-        else
-          do l = 1, nlayer
-            pdqc(ig,l,iq) = (1./masse(l)) * ( zmflux(l)*(zqm(l,iq)-zq(l,iq)) - zmflux(l+1)*(zqm(l+1,iq)-zq(l,iq))&
-                            + zcondicea(ig,l)*zq(l,iq) )
-          end do
-        end if
-      end do
-    end if ! if
-  end do  ! loop on ig
-!----------------------------------------------------------------------------------------------------------------------!
-! 6. CO2 snow / clouds scheme
-!----------------------------------------------------------------------------------------------------------------------!
-  call co2snow(ngrid, nlayer, ptimestep, emisref, condsub, pplev, zcondicea, zcondices, zdqssed_co2*ptimestep, &
-               pemisurf)
-!----------------------------------------------------------------------------------------------------------------------!
-! 7. Extra special case for surface temperature tendency pdtsrfc:
-!      We want to fix the south pole temperature to CO2 condensation temperature.
-!----------------------------------------------------------------------------------------------------------------------!
-#ifndef MESOSCALE 
-  if (caps.and.(obliquit<27.)) then
-    ! check if last grid point is the south pole
-    if (abs(latitude(ngrid)-(-pi/2.))<1.e-5) then
-!     NB: Updated surface pressure, at grid point 'ngrid', is ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
-!     write(*,*)"co2condens: South pole: latitude(ngrid)=", latitude(ngrid)
-      ztcondsol(ngrid) = 1./(bcond-acond*log(.01*vmr_co2(ngrid,1) * (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
-      pdtsrfc(ngrid) = (ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
-    end if
-  end if
-#endif
-!======================================================================================================================!
-! END =================================================================================================================!
-!======================================================================================================================!
-  end subroutine co2condens4micro
-
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-
-!======================================================================================================================!
-! SUBROUTINE: Van-Leer scheme =========================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!  Operateur de moyenne inter-couche pour calcul de transport type Van-Leer " pseudo amont " dans la verticale
-!----------------------------------------------------------------------------------------------------------------------!
-! Comments:
-!----------
-!  q,w are input arguments for the s-pg ....
-!----------------------------------------------------------------------------------------------------------------------!
-! Paper:
-!-------
-!  Van-Leer (1977), "Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection"
-!----------------------------------------------------------------------------------------------------------------------!
-  subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
-
-  implicit none
-!----------------------------------------------------------------------------------------------------------------------!
-! VARIABLE DECLARATION
-!----------------------------------------------------------------------------------------------------------------------!
-! Input arguments:
-!-----------------
-  integer, intent(in) :: &
-     nlayer ! number of layers
-
-  real, intent(in) :: &
-     pente_max,     &! coefficient, pente_max = 2 advised
-     masse(nlayer), &! masse layer Dp/g (kg)
-     q(nlayer)       ! quantity of tracer
-!----------------------------------------------------------------------------------------------------------------------!
-! In-Output arguments:
-!---------------------
-  real, intent(inout) :: &
-     w(nlayer+1) ! masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
-!----------------------------------------------------------------------------------------------------------------------!
-! Output arguments:
-!------------------
-  real, intent(out) :: &
-     qm(nlayer+1) ! quantity of tracer after Van-Leer scheme
-!----------------------------------------------------------------------------------------------------------------------!
-! Locals variables:
-!------------------
-  integer :: &
-     l, &! loop on nlayer
-     m   ! index
-
-  real :: &
-     dzqmax,      &! maximum of dzq between two adjacent layers
-     sigw,        &!
-     Mtot,        &!
-     MQtot,       &!
-     dzq(nlayer), &!
-     dzqw(nlayer),&!
-     adzqw(nlayer) !
-!======================================================================================================================!
-! BEGIN ===============================================================================================================!
-!======================================================================================================================!
-! 1. On oriente tout dans le sens de la pression: w > 0 when down
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 2, nlayer
-    dzqw(l) = q(l-1) - q(l)
-    adzqw(l) = abs(dzqw(l))
-  end do
-
-  do l = 2, nlayer-1
-    if(dzqw(l)*dzqw(l+1)>0.) then
-      dzq(l) = 0.5 * (dzqw(l)+dzqw(l+1))
-    else
-      dzq(l) = 0.
-    end if
-
-    dzqmax = pente_max * min(adzqw(l), adzqw(l+1))
-
-    dzq(l) = sign(min(abs(dzq(l)),dzqmax), dzq(l))
-  end do
-
-  dzq(1)=0.
-  dzq(nlayer)=0.
-
-  do l = 1, nlayer-1
-!----------------------------------------------------------------------------------------------------------------------!
-! 2.1. Regular scheme (transfered mass < layer mass)
-!----------------------------------------------------------------------------------------------------------------------!
-    if (w(l+1)>0. .and. w(l+1)<=masse(l+1)) then
-      sigw = w(l+1) / masse(l+1)
-      qm(l+1) = (q(l+1) + 0.5*(1.-sigw)*dzq(l+1))
-    else if (w(l+1)<=0. .and. -w(l+1)<=masse(l)) then
-      sigw = w(l+1) / masse(l)
-      qm(l+1) = (q(l) - 0.5*(1.+sigw)*dzq(l))
-!----------------------------------------------------------------------------------------------------------------------!
-! 2.2. Extended scheme (transfered mass > layer mass)
-!----------------------------------------------------------------------------------------------------------------------!
-    else if (w(l+1)>0.) then
-      m = l+1
-      Mtot = masse(m)
-      MQtot = masse(m)*q(m)
-
-      do while ((m<nlayer).and.(w(l+1)>(Mtot+masse(m+1))))
-        m = m+1
-        Mtot = Mtot + masse(m)
-        MQtot = MQtot + masse(m)*q(m)
-      end do
-
-      if (m<nlayer) then
-        sigw = (w(l+1)-Mtot) / masse(m+1)
-        qm(l+1) = (1/w(l+1))*( MQtot + (w(l+1)-Mtot)* (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
-      else
-        w(l+1) = Mtot
-        qm(l+1) = Mqtot / Mtot
-        call abort_physic('co2condens4micro', 'top layer is disapearing !', 1)
-      end if
-!----------------------------------------------------------------------------------------------------------------------!
-    else ! if(w(l+1).lt.0)
-      m = l-1
-      Mtot = masse(m+1)
-      MQtot = masse(m+1)*q(m+1)
-      if (m>0) then ! because some compilers will have problems evaluating masse(0)
-        do while ((m>0).and.(-w(l+1)>(Mtot+masse(m))))
-          m = m-1
-          Mtot = Mtot + masse(m+1)
-          MQtot = MQtot + masse(m+1)*q(m+1)
-          if (m==0) then
-            call abort_physic('co2condens_mod4micro','vl1d',1)
-          end if
-        end do
-      end if
-
-      if (m>0) then
-        sigw = (w(l+1)+Mtot) / masse(m)
-        qm(l+1) = (-1/w(l+1)) * ( MQtot + (-w(l+1)-Mtot) * (q(m)-0.5*(1.+sigw)*dzq(m)) )
-      else
-        qm(l+1) = (-1/w(l+1)) * (MQtot + (-w(l+1)-Mtot)*qm(1))
-      end if
-    end if
-  end do ! l = 1, nlayer-1
-!======================================================================================================================!
-! END =================================================================================================================!
-!======================================================================================================================!
-  end subroutine vl1d
-
-end module co2condens_mod4micro
Index: trunk/LMDZ.MARS/libf/phymars/co2sat.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2sat.F	(revision 2348)
+++ trunk/LMDZ.MARS/libf/phymars/co2sat.F	(revision 2349)
@@ -1,41 +1,51 @@
+      SUBROUTINE co2sat(naersize,t,p,psat)
+c       SUBROUTINE co2sat(naersize,t,p,qsat) JA
+      IMPLICIT NONE
+
 c=======================================================================
-c     SUBROUTINE co2sat
-c-----------------------------------------------------------------------
-c     Aim:
-c     ----
-c     Compute saturated steam pressure (from James et al, 1992)
+c
+c
+c  now:  straight psat of CO2 (or qsat of CO2 but need of mmean)
+c
 c=======================================================================
-      subroutine co2sat(naersize, t, psat)
-     
-      implicit none
-c-----------------------------------------------------------------------
-c     VARIABLES
-c-----------------------------------------------------------------------
-c     Inputs:
-c     ------- 
-      integer, intent(in) ::
-     &   naersize ! dimension of tables t and psat
-     
-      real, intent(in) ::
-     &   t(naersize) ! temperature table
 
-c     Output:
-c     -------
-      real, intent(out) ::
-     &   psat(naersize) ! Saturated steam pressure (Pa)
+c   declarations:
+c   -------------
+c   arguments:
+c   ----------
 
-c     Local:
-c     ------
-      integer ::
-     &   i ! loop on naersize
-c=======================================================================
-c===== BEGIN
-c=======================================================================
-      do i = 1, naersize
-        psat(i) = 1.382 * 1e12 * exp(-3182.48/t(i))
-      end do
-c=======================================================================
-c===== END
-c=======================================================================
+c   INPUT
+      integer naersize 
+      real t(naersize) , p(naersize)
+c   OUTPUT
+c      real qsat(naersize) JA
+      real psat(naersize)
+
+c   local:
+c   ------
+      INTEGER i
+      REAL r2,r3,r4 , To, es
+      SAVE r2,r3,r4
+      DATA r2,r3,r4/611.14,21.875,7.66/
+      SAVE To
+      DATA To/273.16/
+          
+      do i=1,naersize
+
+
+c        pression de vapeur saturante (James et al. 1992):
+
+          psat(i)  = 1.382 * 1e12 * exp(-3182.48/t(i)) !; (Pa)
+
+c         OR:
+
+c         qsat(i) = psat/p(i)*44.01/mmean ! Need of updated information on mmean
+c         qsat(i) = max(qsat(i), 1.e-30)
+
+
+      enddo
+c      qsat=psat JA
+          
+
       RETURN
       END
Index: trunk/LMDZ.MARS/libf/phymars/co2snow.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 2348)
+++ trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 2349)
@@ -1,141 +1,138 @@
-c=======================================================================
-c     Program for simulate the impact of the CO2 snow fall on
-c     the surface infrared emission (emissivity)
-c-----------------------------------------------------------------------
-c     Algorithme:
-c     1 - Initialization
-c     2 - Compute the surface emissivity
-c-----------------------------------------------------------------------
-c     Author: F. Forget 
-c-----------------------------------------------------------------------
-c     Reference paper:
-c     Francois Forget and James B. Pollack
-c
-c     "Thermal infrared observations of the condensing Martian polar
-c     caps: CO2 ice temperatures and radiative budget"
-c
-c     JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. E7, PAGES 16,
-c     865-16,879, JULY 25, 1996
-c=======================================================================
+      SUBROUTINE co2snow (ngrid,nlayer,ptimestep,emisref,condsub
+     &         ,pplev,pcondicea,pcondices,pfallice,pemisurf)
 
-      SUBROUTINE co2snow(ngrid, nlayer, ptimestep, emisref, condsub,
-     &                   pplev, pcondicea, pcondices, pfallice,
-     &                   pemisurf)
-      
       use surfdat_h, only: iceradius, dtemisice
       use geometry_mod, only: latitude ! grid point latitudes (rad)
       use time_phylmdz_mod, only: daysec
+      IMPLICIT NONE
 
-      IMPLICIT NONE
+c=======================================================================
+c     Program for simulate the impact of the CO2 snow fall on
+c     the surface infrared emission (emissivity)  and on 
+c     the airborne dust
+c     F.Forget 1996
+c=======================================================================
+
+c Declarations
+c ------------
 
 #include "callkeys.h"
 
+c     Arguments
+c     ---------
 
-c=======================================================================
-c     Variables
-c=======================================================================
-c     Inputs:
-c     -------
-      integer, intent(in) ::
-     &   ngrid, ! number of atmospheric columns
-     &   nlayer ! number of atmospheric layers
-                        
-      real, intent(in) ::
-     &   ptimestep,               ! timestep of the physics (s)
-     &   emisref(ngrid),          ! grd or ice emissivity without snow 
-     &   pplev(ngrid,nlayer+1),   ! interlayer pressure (Pa)
-     &   pcondicea(ngrid,nlayer), ! CO2 condensation rate in layers
-                                  !   (kg/m2/s)
-     &   pcondices(ngrid),        ! CO2 condensation rate on the surface
-                                  !   (kg/m2/s)        
-     &   pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
-                           
-      logical, intent(in) ::
-     &   condsub(ngrid) ! true if there is CO2 condensation or
-                        ! sublimation in the column
+      INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      REAL,INTENT(IN) :: ptimestep ! timestep of the physics (s)
+      REAL,INTENT(IN) :: emisref(ngrid) ! grd or ice  emissivity without snow
+      logical,intent(in) :: condsub(ngrid) ! true if there is CO2 condensation
+                                           ! or sublimation in the column
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! interlayer pressure (Pa)
+      REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! CO2 condensation rate (kg/m2/s)
+      REAL,INTENT(IN) :: pcondices(ngrid) ! CO2 condensation rate (kg/m2/s)
+                                          ! on the surface
+      REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
 
-c     Output:
-c     -------
-      real, intent(out) ::
-     &   pemisurf(ngrid) ! surface emissivity
+      REAL,INTENT(OUT) :: pemisurf(ngrid) ! surface emissivity
 
-c     local:
-c     ------
-      integer ::
-     &   l,
-     &   ig,
-     &   icap
+c     local
+c     -----
+      integer ig , l , icap
 
-      real ::
-     &   sumdaer,
-     &   zdemisurf
+      REAL zdemisurf ,dtemis
+      REAL sumdaer
 
-      real, parameter ::
-     &   alpha = 0.45
+c     saved
+c     -----
+      REAL,save :: Kscat(2), scaveng
+      LOGICAL,SAVE :: firstcall=.true.
 
-c     saved:
-c     ------
-      real, save ::
-     &   Kscat(2) ! Kscat: coefficient for decreasing the surface 
-c                 !        emissivity
-c                 ! 
-c                 ! Kscat = (0.001/3.) * alpha / iceradius
-c                 !          with 0.3 < alpha < 0.6
-c                 !
-c                 ! alpha set to 0.45 (coeff from emis = f (tau)) and
-c                 ! iceradius the mean radius of the scaterring
-c                 ! particles (200.e-6 < iceradius < 10.e-6)
-c                 !
-c                 ! (2) = N and S hemispheres
+c --------------
+c Initialisation
+c --------------
 
-      logical, save ::
-     &   firstcall = .true.
-c=======================================================================
-c BEGIN
-c=======================================================================
-c 1 - Initialization
-c=======================================================================
+      ! AS: firstcall OK absolute
       if (firstcall) then
-        Kscat(1) = (0.001/3.) * alpha / iceradius(1)
-        Kscat(2) = (0.001/3.) * alpha / iceradius(2)
-                   
-        firstcall = .false.
+
+c   Kscat : coefficient for decreasing the surface emissivity
+c           =(0.001/3.)*alpha/iceradius ,
+c           with 0.3< alpha < 0.6, set to 0.45 (coeff from emis = f (tau))
+c           and iceradius the mean radius of the
+c           scaterring particles  (200.e-6<iceradius<10.e-6)
+
+            Kscat(1)=(0.001/3.)*0.45/iceradius(1)
+            Kscat(2)=(0.001/3.)*0.45/iceradius(2)
+
+c          Scavenging Ratio (dust concentration in the air / in the snow)
+           scaveng = 100.0
+           
+c          Collision Scavenging coefficient  (m2.kg-1)
+c          Csca  = 2.3  ! not used yet !!!!!!!!!!!
+           firstcall = .false.
+
       end if
-c=======================================================================
-c 2 - Compute the surface emissivity
-c=======================================================================
-      do ig = 1, ngrid
-        if (condsub(ig)) then
-          if (latitude(ig).lt.0.) then
-              icap = 2 ! Southern hemisphere
-          else
-              icap = 1 ! Northern Hemisphere
-          end if
 
-          ! compute zdemisurf using an integrated form for numerical
-          ! safety instead
-          zdemisurf =
-     &     (emisref(ig) - pemisurf(ig)) / (dtemisice(icap) * daysec)
-     &     +
-     &     ( emisref(ig) * 
-     &              ( (pemisurf(ig)/emisref(ig))**(-3) +
-     &                3.*Kscat(icap)*pfallice(ig,1)*ptimestep )**(-1/3.)
-     &       - pemisurf(ig)
-     &      ) / ptimestep
-          
-          pemisurf(ig) = pemisurf(ig) + zdemisurf * ptimestep 
-          if (pemisurf(ig).lt.0.1) then
-            write(*,*)'ds co2snow: emis < 0.1 !!!'
-            write(*,*)'ig =' , ig
-            write(*,*)'pemisurf(ig)', pemisurf(ig)
-            write(*,*)'zdemisurf*ptimestep', zdemisurf*ptimestep
-          end if
-        else ! if condsub(ig) is false
+
+c     LOOP on grid points
+c     -------------------
+      do ig=1,ngrid 
+         if (condsub(ig)) then
+
+c         IF (scavenging) then
+c          Airborne Dust 
+c          -------------
+c          sumdaer=0.
+c          do l=nlayer, 1, -1
+c             pdaerosol(ig,l)= -paerosol(ig,l,1)*
+c    &              (1-exp(-scaveng*pcondicea(ig,l)*ptimestep*g/
+c    &               (pplev(ig,l)-pplev(ig,l+1))))/ptimestep  
+
+c    &           - Csca*paerosol(ig,l,1) ! Scavenging by collision
+c    &           * 0.5*(pfallice(ig,l)) ! not included
+
+c             test to avoid releasing to much dust when subliming:
+c             if(pdaerosol(ig,l).gt.-sumdaer)pdaerosol(ig,l)=-sumdaer 
+c             sumdaer=sumdaer + pdaerosol(ig,l)
+
+c            if (-pdaerosol(ig,l)*ptimestep.gt.paerosol(ig,l,1)) then
+c                write(*,*) 'ds co2snow: aerosol < 0.0 !!!'
+c                write(*,*) 'ig =' , ig
+c            end if
+c          end do
+c         END IF
+
+c          Surface emissivity
+c          ------------------
+c   dtemis: Time scale for increasing the ice emissivity
+
+           IF(latitude(ig).LT. 0.) THEN
+              icap=2 ! Southern hemisphere
+           ELSE
+              icap=1 ! Northern Hemisphere
+           ENDIF
+
+           zdemisurf = 
+     &    (emisref(ig)-pemisurf(ig))/(dtemisice(icap)*daysec) 
+c Using directly the diferential equation:
+c    &       -Kscat(icap)*emisref(ig)*
+c    &      (pemisurf(ig)/emisref(ig))**4 *pfallice(ig,1) 
+c Using an integrated form for numerical safety instead
+     & +(emisref(ig)* ((pemisurf(ig)/emisref(ig))**(-3)+3.*Kscat(icap)*
+     & pfallice(ig,1)*ptimestep)**(-1/3.) -pemisurf(ig))/ptimestep   
+
+
+           pemisurf(ig) = pemisurf(ig) + zdemisurf*ptimestep 
+
+           if (pemisurf(ig).lt.0.1) then
+                 write(*,*) 'ds co2snow: emis < 0.1 !!!'
+                 write(*,*) 'ig =' , ig
+                 write(*,*)'pemisurf(ig)',pemisurf(ig)
+                 write(*,*) 'zdemisurf*ptimestep',zdemisurf*ptimestep
+           end if
+         else
            pemisurf(ig) = emisref(ig)
-        end if
-      end do ! ngrid
-c=======================================================================
-c END
-c=======================================================================
+         end if
+      end do
+
       return
       end 
Index: trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F	(revision 2349)
+++ trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F	(revision 2349)
@@ -0,0 +1,671 @@
+      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 ( (Ic_rice.ne.Ic_rice) ! will be true if it is Nan  
+     &               .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.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90	(revision 2348)
+++ 	(revision )
@@ -1,624 +1,0 @@
-!======================================================================================================================!
-! Module: Scheme of co2 cloud formation ===============================================================================!
-!----------------------------------------------------------------------------------------------------------------------!
-! Authors: Joaquim Audouard, Constantino Listowski, Anni Määttänen
-! Date: 09/2016
-!----------------------------------------------------------------------------------------------------------------------!
-! Contains subroutines:
-!     - improvedco2clouds_mod: nucleation
-!
-!     - density_co2_ice: compute density of co2 ice particle
-!======================================================================================================================!
-module improvedco2clouds_mod
-
-implicit none
-
-contains
-!======================================================================================================================!
-! SUBROUTINE: improvedco2clouds =======================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!  This routine is used to form CO2 clouds when a parcel of the GCM is saturated.
-!----------------------------------------------------------------------------------------------------------------------!
-! Comments:
-!----------
-!  Adaptation for CO2 clouds based on improvedclouds_mod.F
-!
-!  It includes the ability to have supersaturation, a computation of the nucleation rates, growthrates and the
-!  scavenging of dust particles by clouds. It is worth noting that the amount of dust is computed using the dust
-!  optical depth computed in aeropacity.F.
-!  That's why the variable called "tauscaling" is used to convert pq(dust_mass) and pq(dust_number), which are relative
-!  quantities, to absolute and realistic quantities stored in zq. This has to be done to convert the inputs into
-!  absolute values, but also to convert the outputs back into relative values which are then used by the sedimentation
-!  and advection schemes.
-!  CO2 ice particles can nucleate on both dust and water ice particles. When CO2 ice is deposited onto a water ice
-!  particles, the particle is removed from the water tracers. Memory of the origin of the co2 particles is kept and
-!  thus the water cycle shouldn't be modified by this.
-!  There is an energy limit to how much co2 can sublimate/condensate. It is defined by the difference of the GCM
-!  temperature with the co2 condensation temperature.
-!
-!                                                /!\ WARNING /!\
-!  No sedimentation of the water ice origin is performed in the microphysical timestep in co2cloud.F.
-!
-!  If meteoritic particles are activated and turn into co2 ice particles, then they will be reversed in the dust
-!  tracers if the cloud sublimates.
-!----------------------------------------------------------------------------------------------------------------------!
-! Paper:
-!-------
-!    see co2cloud.F90
-!----------------------------------------------------------------------------------------------------------------------!
-! Algorithm:
-!-----------
-!   0. Firstcall
-!     0.1. Bonus: meteoritic component, extract data
-!   1. Initialization
-!   2. Compute saturation
-!   3. Bonus: additional meteoritic particles for nucleation
-!   4. Actual microphysics: Main loop over the GCM's grid
-!     4.1 Nucleation
-!     4.2. Ice growth: scheme for radius evolution
-!     4.3 Dust cores releasing if no more co2 ice
-!   5. Get cloud tendencies
-!======================================================================================================================!
-  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, rb_cldco2, sigma_iceco2, dev2)
-
-  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
-
-  include "callkeys.h"
-  include "microphys.h"
-!----------------------------------------------------------------------------------------------------------------------!
-! VARIABLES DECLARATION
-!----------------------------------------------------------------------------------------------------------------------!
-! Input arguments:
-!-----------------
-  integer, intent(in) :: &
-     nq,    &! number of tracers
-     ngrid, &! number of point grid
-     nlay    ! number of layer
-
-  real, intent(in) :: &
-     microtimestep,           &! physics time step (s)
-     pplay(ngrid,nlay),       &! mid-layer pressure (Pa)
-     pplev(ngrid,nlay+1),     &! inter-layer pressure (Pa)
-     pteff(ngrid,nlay),       &! temperature at the middle of the layers (K)
-     sum_subpdt(ngrid,nlay),  &! tendency on temperature from previous physical parametrizations
-     pqeff(ngrid,nlay,nq),    &! tracers (kg/kg)
-     tauscaling(ngrid),       &! convertion factor for qdust and Ndust
-     sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers before condensation (kg/kg.s-1)
-
-  real,  intent(in) :: &
-     sigma_iceco2 ! Variance of the co2 ice and CCN distributions
-
-  double precision, intent(in) :: &
-     rb_cldco2(nbinco2_cld+1), & ! boundary values of each rad_cldco2
-     dev2
-!----------------------------------------------------------------------------------------------------------------------!
-! Output arguments:
-!------------------
-  real, intent(out) :: &
-     subpdtcloudco2(ngrid,nlay),  &! tendency on tracers due to CO2 condensation (K/s)
-     subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers due to CO2 condensation (kg/kg.s-1)
-!----------------------------------------------------------------------------------------------------------------------!
-! Local:
-!-------
-!----1) Parameters:
-!------------------
-  integer, parameter :: &
-! ---Meteoritic flux input file
-     nbin_meteor = 100, &! Dimension 1
-     nlev_meteor = 130, &! Dimension 2
-     uMeteor = 666,     &! File unit
-! ---Latent heat computation
-     l0 = 595594d0,     &! coeff from: Azreg-Aïnou (2005)
-     l1 = 903.111d0,    &!   Title: "Low-temperature data for carbon dioxide"
-     l2 = -11.5959d0,   &!   Pulication: eprint arXiv:1403.4403
-     l3 = 0.0528288d0,  &!
-     l4 = -0.000103183d0 !
-
-  real, parameter :: &
-     threshold = 1e-30 ! limit value
-
-  character(len=20), parameter:: &
-     file_meteoritic_flux = 'Meteo_flux_Plane.dat'
-!----------------------------------------------------------------------------------------------------------------------!
-!----2) Saved:
-!-------------
-  real, save :: &
-     sigma_ice  ! Variance of the h2o ice and CCN distributions
-
-  double precision, save :: &
-     meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor
-     dev3                              ! 1. / ( sqrt(2.) * sigma_ice )
-
-  logical, save :: &
-     firstcall = .true. ! Used to compute saved variables
-!----------------------------------------------------------------------------------------------------------------------!
-!----3) Variables:
-!-----------------
-  integer :: &
-     ig,     &! loop on ngrid
-     l,      &! loop on nlay
-     i,      &! loop on nbinco2
-! ---Variables for meteoritic flux
-     ibin,   &! loop on nbin_meteor
-     nelem,  &! nb of points during interpolation of meteoritic flux
-     lebon1, &! index where P_meteor is the nearest of pplev(ig,l)
-     lebon2, &! index where P_meteor is the nearest of pplev(ig,l+1)
-     read_ok  ! file uMeteor iostat
-
-  real :: &
-    masse(ngrid,nlay),     &! mass layer (kg.m-2)
-    rice(ngrid,nlay),      &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
-    zq(ngrid,nlay,nq),     &! local value of tracers (kg/kg)
-    zq0(ngrid,nlay,nq),    &! local init value of tracers (kg/kg)
-    zt(ngrid,nlay),        &! local value of temperature (K)
-    zqsat(ngrid,nlay),     &! saturation vapor pressure for CO2 (K)
-    tcond(ngrid,nlay),     &! condensation temperature of CO2 (K)
-    lw,                    &! Latent heat of sublimation (J.kg-1)
-    rdust(ngrid,nlay),     &! Dust geometric mean radius (m)
-    rhocloud(ngrid,nlay),  &! Cloud density (kg.m-3)
-    rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
-
-  real(kind=8) :: &
-    derf ! Error function
- 
-  double precision :: &
-    dMice,                             &! mass of condensed ice
-    facteurmax,                        &! for energy limit on mass growth
-    pco2,                              &! Co2 vapor partial pressure (Pa)
-    satu,                              &! Co2 vapor saturation ratio over ice
-    Mo,                                &! mass of aerosol particles
-    No,                                &! number of aerosol particles
-    Rn,                                &! logarithm of rdust/rice
-    Rm,                                &! Rn * variance of ice and CCN distribution
-    n_derf,                            &! derf( (rb_cldco2(1)+Rn) *dev3)
-    m_derf,                            &! derf( (rb_cldco2(1)+Rm) *dev2)
-    mem_Mccn_co2(ngrid,nlay),          &! Memory of CCN mass of H2O and dust used by CO2
-    mem_Mh2o_co2(ngrid,nlay),          &! Memory of H2O mass integred into CO2 crystal
-    mem_Nccn_co2(ngrid,nlay),          &! Memory of CCN number of H2O and dust used by CO2
-    n_aer(nbinco2_cld),                &! Radius used by the microphysical scheme (m)
-    m_aer(nbinco2_cld),                &! number concentration V-1 of particle/each size bin
-    n_aer_h2oice(nbinco2_cld),         &! mass mixing ratio of particle/each size bin
-    m_aer_h2oice(nbinco2_cld),         &! Same - for CO2 nucleation
-    rad_h2oice(nbinco2_cld),           &! Same - for CO2 nucleation
-    Ic_rice,                           &! Mass transfer rate CO2 ice crystal
-    ratioh2o_ccn,                      &! 1./(zq(ig,l,igcm_h2o_ice)  + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
-    vo2co2,                            &! volume of co2 ice particle
-    dN,                                &! number of particle of dust used as ccn
-    dM,                                &! mass of dN
-    dNh2o,                             &! number of particle of h2o ice used as ccn
-    dMh2o,                             &! mass of dNh2o
-    dNN,                               &! min(dN,zq(ig,l,igcm_dust_number))
-    dMM,                               &! min(dM,zq(ig,l,igcm_dust_mass))
-    dNNh2o,                            &! min(dNNh2o,zq(ig,l,igcm_ccn_number))
-    dMh2o_ice,                         &! min(dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn, zq(ig,l,igcm_h2o_ice))
-    dMh2o_ccn,                         &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
-    rate(nbinco2_cld),                 &! nucleation rate
-    rateh2o(nbinco2_cld),              &! nucleation rate for h2o
-    rho_ice_co2T(ngrid,nlay),          &! density of co2 ice
-    riceco2(ngrid,nlay),               &! CO2 ice mean radius (m)
-    vrat_cld,                          &! Volume ratio
-    Proba,                             &! 1.0 - dexp(-1.*microtimestep*rate(i))
-    Probah2o,                          &! 1.0 - dexp(-1.*microtimestep*rateh2o(i))
-    mtemp(nbinco2_cld),                &! sum(meteor(lebon1:lebon2,ibin))
-    pression_meteor(nlev_meteor),      &! pressure from meteoritic flux input file
-    ltemp1(nlev_meteor),               &! abs(pression_meteor(:)-pplev(ig,l))
-    ltemp2(nlev_meteor),               &! abs(pression_meteor(:)-pplev(ig,l+1))
-    meteor_ccn(ngrid,nlay,nbinco2_cld)  ! nbinco2_cld = 100
-
-  logical :: &
-    file_ok ! test if meteoritic input file exists
-!======================================================================================================================!
-! BEGIN ===============================================================================================================!
-!======================================================================================================================!
-! 0. Firstcall
-!----------------------------------------------------------------------------------------------------------------------!
-  if (firstcall) then
-    firstcall = .false.
-
-!   Volume of a co2 molecule (m3)
-    vo1co2 = m0co2 / dble(rho_ice_co2)
-
-!   Variance of the ice and CCN distributions
-    sigma_ice = sqrt(log(1.+nuice_sed))
-    dev3 = 1. / ( sqrt(2.) * sigma_ice )
-!----------------------------------------------------------------------------------------------------------------------!
-! 0.1. Bonus: meteoritic component, extract data
-!----------------------------------------------------------------------------------------------------------------------!
-    if (meteo_flux) then
-      ! Check if file exists
-      inquire(file=trim(datadir)//'/'//file_meteoritic_flux, exist=file_ok)
-      if (.not. file_ok) then
-        call abort_physic("CO2clouds", 'file '//file_meteoritic_flux//' should be in'//trim(datadir), 1)
-      end if
-
-      ! open file
-      open(unit=uMeteor,file=trim(datadir)//'/'//file_meteoritic_flux, FORM='formatted')
-
-      !skip 1 line
-      read(uMeteor,*)
-
-      ! extract pressure_meteor
-      do i = 1, nlev_meteor
-        read(uMeteor,*)pression_meteor(i)
-      end do
-
-      !skip 1 line
-      read(uMeteor,*)
-
-      ! extract meteor flux
-      do i = 1, nlev_meteor
-        ! les mêmes 100 bins size que la distri nuclea : on touche pas
-        do ibin = 1, nbin_meteor
-          read(uMeteor,'(F12.6)') meteor(i,ibin)
-        end do
-      end do
-
-      ! close file
-      close(uMeteor)
-    end if ! of if meteo_flux
-  end if ! firstcall
-!----------------------------------------------------------------------------------------------------------------------!
-! 1. Initialization
-!----------------------------------------------------------------------------------------------------------------------!
-  rdust(:,:) = 0.
-  meteor_ccn(:,:,:) = 0.
-  rice(:,:) = 1.e-8
-  riceco2(:,:) = 1.e-11
-
-  ! Initialize the tendencies
-  subpdqcloudco2(:,:,:) = 0.
-  subpdtcloudco2(:,:) = 0.
-
-  ! pteff temperature layer; sum_subpdt dT.s-1
-  zt(:,:) = pteff(:,:) + sum_subpdt(:,:) * microtimestep
-
-  ! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
-  zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep
-
-  where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold
-
-  zq0(:,:,:) = zq(:,:,:)
-!----------------------------------------------------------------------------------------------------------------------!
-! 2. Compute saturation
-!----------------------------------------------------------------------------------------------------------------------!
-  call co2sat(ngrid*nlay,zt,zqsat)
-  call tcondco2(ngrid,nlay,pplay, zq(:,:,igcm_co2) + zq(:,:,igcm_co2_ice),tcond)
-!----------------------------------------------------------------------------------------------------------------------!
-! 3. Bonus: additional meteoritic particles for nucleation
-!----------------------------------------------------------------------------------------------------------------------!
-  ! TODO: instead of intepolation, used only the nearest pplev(ig,l)
-  if (meteo_flux) then
-    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
-
-        do ibin = 1, nbin_meteor
-          mtemp(ibin) = sum(meteor(lebon1:lebon2,ibin))
-        end do
-
-        ! Par kg air csi par m carre, x epaisseur/masse pour par kg/air. Check original unit with J. Plane
-        meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l)
-      end do
-    end do
-  end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 4. Actual microphysics: Main loop over the GCM's grid
-!----------------------------------------------------------------------------------------------------------------------!
-  do l = 1, nlay
-    do ig = 1, ngrid
-
-      ! Get the partial pressure of co2 vapor and its saturation ratio
-      pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) ! mco2 (kg/mol) => mmean (g/mol)
-      satu = pco2 / zqsat(ig,l)
-
-      !T-dependant CO2 ice density
-      call density_co2_ice(zt(ig,l), rho_ice_co2T(ig,l))
-
-      vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
-      rho_ice_co2 = rho_ice_co2T(ig,l)
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.1 Nucleation
-!----------------------------------------------------------------------------------------------------------------------!
-      ! if there is condensation
-      if ( satu >= 1 ) then
-        call updaterccnCO2(zq(ig,l,igcm_dust_mass), zq(ig,l,igcm_dust_number), rdust(ig,l), tauscaling(ig))
-
-        ! Expand the dust moments into a binned distribution
-        n_aer(:) = 0d0 ! number of aerosol particles
-        m_aer(:) = 0d0 ! mass of aerosol particles
-
-        No = zq(ig,l,igcm_dust_number) * tauscaling(ig)
-
-        Mo = (4./3.) * pi * rho_dust * No * rdust(ig,l)**3 *dexp(9.*nuiceco2_ref/2.)
-
-        if (No > threshold) then
-          Rn = -dlog(rdust(ig,l))
-
-          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
-            m_aer(i) = -0.5 * Mo * m_derf
-
-            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
-          end do
-
-          ! 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)**3
-             end do
-          end if
-
-          ! 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) + threshold
-
-            ! Total mass of H20 crystals,CCN included
-            No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
-
-            Rn = -dlog(rice(ig,l))
-
-            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)
-            end do
-          end if
-
-          ! 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 - dexp(-1.*microtimestep*rate(i))
-            dN = dN + n_aer(i) * Proba
-            dM = dM + m_aer(i) * Proba
-          end do
-
-          if (co2useh2o) then
-            do i = 1, nbinco2_cld
-              Probah2o = 1.0 - dexp(-1.*microtimestep*rateh2o(i))
-              dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
-              dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
-            end do
-          end if
-
-          ! Now increment CCN tracers and update dust tracers
-          dNN = min(dN,zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
-          dMM = min(dM,zq(ig,l,igcm_dust_mass))  ! idem dans le min
-
-          zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
-
-          zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
-
-          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
-
-          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
-
-          ! 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_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 = dMh2o * zq(ig,l,igcm_h2o_ice) * ratioh2o_ccn
-            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
-          end if ! of if co2useh2o
-        end if ! of if No > 1e-30
-      end if ! of is satu > 1
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.2. Ice growth: scheme for radius evolution
-!----------------------------------------------------------------------------------------------------------------------!
-!    We trigger crystal growth if and only if there is at least one nuclei (N>1). Indeed, if we are supersaturated
-!      and still don't have at least one nuclei, we should better wait to avoid unrealistic value for nuclei radius
-!      and so on for cases that remain negligible.
-!----------------------------------------------------------------------------------------------------------------------!
-      ! we trigger crystal growth
-      if (zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) + threshold >= 1.) then
-
-        call updaterice_microco2(zq(ig,l,igcm_co2_ice), zq(ig,l,igcm_ccnco2_mass), zq(ig,l,igcm_ccnco2_number), &
-                                 tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
-
-        Ic_rice = 0.
-
-        ! J.kg-1
-        lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4
-
-        facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw)
-
-        ! 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)
-
-        ! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance !
-        if (isnan(Ic_rice) .or. Ic_rice == 0.) then
-          Ic_rice = 0.
-          subpdtcloudco2(ig,l) = -sum_subpdt(ig,l)
-          dMice = 0
-        else
-          ! Kg par kg d'air, >0 si croissance !
-          ! kg.s-1 par particule * nb particule par kg air*s = kg par kg air
-          dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig)
-
-          ! 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
-          dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
-          dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
-
-          ! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s
-          subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep
-
-          !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
-        end if
-!----------------------------------------------------------------------------------------------------------------------!
-! 4.3 Dust cores releasing if no more co2 ice
-!----------------------------------------------------------------------------------------------------------------------!
-        ! On sublime tout
-        if (zq(ig,l,igcm_co2_ice) <= threshold) then
-
-          if (co2useh2o) then
-
-            if (mem_Mccn_co2(ig,l) > 0) then
-              zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + mem_Mccn_co2(ig,l)
-            end if
-
-            if (mem_Mh2o_co2(ig,l) > 0) then
-              zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + mem_Mh2o_co2(ig,l)
-            end if
-
-            if (mem_Nccn_co2(ig,l) > 0) then
-              zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + mem_Nccn_co2(ig,l)
-            end if
-
-          end if
-
-          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.
-        end if !of if co2_ice < threshold
-      end if ! of if NCCN > 1
-    end do ! of ig loop
-  end do ! of nlayer loop
-!----------------------------------------------------------------------------------------------------------------------!
-! 5. Get cloud tendencies
-!----------------------------------------------------------------------------------------------------------------------!
-  subpdqcloudco2(:,:,igcm_co2) =   ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep
-
-  subpdqcloudco2(:,:,igcm_co2_ice) =  ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep
-
-  subpdqcloudco2(:,:,igcm_ccnco2_mass) =  ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep
-
-  subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep
-
-  subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) )  / microtimestep
-
-  subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) )   / microtimestep
-
-  if (co2useh2o) then
-    subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) )  / microtimestep
-
-    subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) )    / microtimestep
-
-    subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
-  end if
-!======================================================================================================================!
-! END =================================================================================================================!
-!======================================================================================================================!
-  end subroutine improvedCO2clouds
-
-
-!**********************************************************************************************************************!
-!**********************************************************************************************************************!
-
-
-!======================================================================================================================!
-! SUBROUTINE: density_co2_ice =========================================================================================!
-!======================================================================================================================!
-! Subject:
-!---------
-!   Compute co2 ice particles density
-!======================================================================================================================!
-  subroutine density_co2_ice(temperature, density)
-
-  implicit none
-
-  double precision, intent(in) :: &
-     temperature
-
-  double precision, intent(out) :: &
-     density
-
-  density = 1000. * (1.72391 - 2.53e-4*temperature - 2.87e-6*temperature*temperature)
-
-  end subroutine density_co2_ice
-
-end module improvedCO2clouds_mod
-
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2348)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2349)
@@ -18,5 +18,4 @@
       use watersat_mod, only: watersat
       use co2condens_mod, only: co2condens
-      use co2condens_mod4micro, only: co2condens4micro
       use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
      &                        mem_Nccn_co2
@@ -252,4 +251,6 @@
 c     Variables used by the CO2 clouds microphysical scheme:
       DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
+      real rsedcloudco2(ngrid,nlayer) !CO2 Cloud sedimentation radius
+      real rhocloudco2(ngrid,nlayer) !co2 Cloud density (kg.m-3)
       real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
       real zcondicea_co2microp(ngrid,nlayer)
@@ -302,5 +303,5 @@
 
 c     Tendancies due to various processes:
-      REAL dqsurf(ngrid,nq)  ! tendency for tracers on surface (Kg/m2/s)
+      REAL dqsurf(ngrid,nq)
       REAL zdtlw(ngrid,nlayer)     ! (K/s)
       REAL zdtsw(ngrid,nlayer)     ! (K/s)
@@ -383,11 +384,9 @@
       REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
       REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
-      REAL mdusttot(ngrid)      ! Total mass of dust tracer (kg/m2)
+      REAL mdusttot(ngrid)    ! Total mass of dust tracer (kg/m2)
       REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
-      REAL mtotco2       ! Total mass of co2, including ice at the surface (kg/m2)
-      REAL vaptotco2     ! Total mass of co2 vapor (kg/m2)
-      REAL icetotco2     ! Total mass of co2 ice (kg/m2)
+      REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
+      REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2) 
       REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
-      REAL NccnCO2tot(ngrid)    ! Total number of ccnCO2 (nbr/m2)
       REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
       REAL rave(ngrid)          ! Mean water ice effective radius (m)
@@ -404,5 +403,5 @@
       REAL DoH_ice(ngrid,nlayer) !D/H ratio
       REAL DoH_surf(ngrid) !D/H ratio surface
-
+                               
       REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
       REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
@@ -484,7 +483,4 @@
       real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
       real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
-C test de conservation de la masse de CO2
-      REAL co2totA
-      REAL co2totB
 
 c entrainment by slope winds above sb-grid scale topography
@@ -500,11 +496,8 @@
 
 c=======================================================================
-c---------- begin added by CM
-      pdq(:,:,:) = 0.
-c---------end added by CM
-
 
 c 1. Initialisation:
 c -----------------
+
 c  1.1   Initialisation only at first call
 c  ---------------------------------------
@@ -563,5 +556,5 @@
            ! starting without startfi.nc and with callsoil
            ! is not yet possible as soildepth default is not defined
-           if (callsoil) then
+           if (callsoil) then                     
               ! default mlayer distribution, following a power law:
               !  mlayer(k)=lay1*alpha**(k-1/2)
@@ -732,5 +725,5 @@
       dsotop(:,:)=0.
       dwatercap(:)=0
-
+      
 #ifdef DUSTSTORM
       pq_tmp(:,:,:)=0
@@ -817,14 +810,4 @@
       enddo
 
-      ! calculates the amount of co2 at the beginning of physics
-      co2totA = 0.
-      do ig=1,ngrid
-        do l=1,nlayer
-           co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
-     &           (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
-     &      +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
-        end do
-        co2totA = co2totA + co2ice(ig)
-      end do
 c-----------------------------------------------------------------------
 c    2. Compute radiative tendencies :
@@ -869,5 +852,5 @@
      &                      pq(1:ngrid,1:nlayer,igcm_o)*
      &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
-
+                 
                  CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
      $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
@@ -919,5 +902,5 @@
 c          ------------------
            ! callradite for the part with clouds
-           clearsky=.false. ! part with clouds for both cases CLFvarying true/false
+           clearsky=.false. ! part with clouds for both cases CLFvarying true/false 
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
      &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
@@ -934,5 +917,5 @@
                ! if
                ! CLFvarying ...) ?? AP ??
-               clearsky=.true. !
+               clearsky=.true. ! 
                CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
      &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
@@ -1136,5 +1119,5 @@
      &                      pdqrds,wspeed,dsodust,dsords,dsotop,
      &                      tauref)
-
+                   
 c      update the tendencies of both dust after vertical transport
          DO l=1,nlayer
@@ -1188,7 +1171,7 @@
            hsummit(:)=14000.
          endif
-         clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
+         clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains 
          nohmons=.false.
-         pdqtop(:,:,:)=0.
+         pdqtop(:,:,:)=0.  
          CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
      &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
@@ -1240,5 +1223,5 @@
         CALL calldrag_noro(ngrid,nlayer,ptimestep,
      &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
-
+ 
         DO l=1,nlayer
           DO ig=1,ngrid
@@ -1326,5 +1309,5 @@
             DO l=1,nlayer
               DO ig=1,ngrid
-                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
               ENDDO
             ENDDO
@@ -1381,5 +1364,5 @@
 
       if(calltherm .and. .not.turb_resolved) then
-
+ 
         call calltherm_interface(ngrid,nlayer,nq,
      $ tracer,igcm_co2,
@@ -1389,5 +1372,5 @@
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
      $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
-
+     
          DO l=1,nlayer
            DO ig=1,ngrid
@@ -1398,8 +1381,8 @@
            ENDDO
         ENDDO
-
+ 
         DO ig=1,ngrid
           q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
-        ENDDO
+        ENDDO      
 
         if (tracer) then
@@ -1510,9 +1493,9 @@
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud,totcloudfrac)
-
+     
 c Temperature variation due to latent heat release
            if (activice) then
-               pdt(1:ngrid,1:nlayer) =
-     &          pdt(1:ngrid,1:nlayer) +
+               pdt(1:ngrid,1:nlayer) = 
+     &          pdt(1:ngrid,1:nlayer) + 
      &          zdtcloud(1:ngrid,1:nlayer)
            endif
@@ -1528,9 +1511,9 @@
            if (hdo) then
 ! increment HDO vapour and ice atmospheric tracers tendencies
-           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
-     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
+           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + 
      &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
-           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
-     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
+           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + 
      &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
            endif !hdo
@@ -1605,13 +1588,15 @@
 c                                nucleation
 c               imicroco2=50     micro-timestep is 1/50 of physical timestep
-        zdqssed_co2(:) = 0.
-
-         IF (co2clouds) THEN
-           call co2cloud(ngrid,nlayer,ptimestep,
+         
+         IF (co2clouds ) THEN
+
+
+            call co2cloud(ngrid,nlayer,ptimestep,
      &           zplev,zplay,pdpsrf,zzlay,pt,pdt,
      &           pq,pdq,zdqcloudco2,zdtcloudco2,
      &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
+     &           rsedcloudco2,rhocloudco2,
      &           rsedcloud,rhocloud,zzlev,zdqssed_co2,
-     &           pdu,pu,zcondicea_co2microp, co2ice)
+     &           pdu,pu,zcondicea_co2microp)
 
 
@@ -1627,23 +1612,24 @@
 ! We need to check that we have Nccn & Ndust > 0
 ! This is due to single precision rounding problems
+              
 ! increment dust tracers tendancies
-               pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
-     &                                 + zdqcloudco2(:,:,igcm_dust_mass)
-
-               pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
-     &                               + zdqcloudco2(:,:,igcm_dust_number)
-
-               pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
-     &                             + zdqcloudco2(:,:,igcm_co2)
-
-               pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
-     &                                 + zdqcloudco2(:,:,igcm_co2_ice)
-
-               pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
-     &                               + zdqcloudco2(:,:,igcm_ccnco2_mass)
-
-               pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
-     &                             + zdqcloudco2(:,:,igcm_ccnco2_number)
-
+               pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)          
+               pdq(1:ngrid,1:nlayer,igcm_dust_number) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
+               pdq(1:ngrid,1:nlayer,igcm_co2) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_co2) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
+               pdq(1:ngrid,1:nlayer,igcm_co2_ice) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_co2_ice) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
+               pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
+               pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
 !Update water ice clouds values as well
              if (co2useh2o) then
@@ -1657,7 +1643,5 @@
      &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
      &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
-
-c     Negative values?
-                where (pq(:,:,igcm_ccn_mass) +
+                where (pq(:,:,igcm_ccn_mass) + 
      &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
                   pdq(:,:,igcm_ccn_mass) = 
@@ -1666,6 +1650,5 @@
      &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
                 end where
-c     Negative values?
-                where (pq(:,:,igcm_ccn_number) +
+                where (pq(:,:,igcm_ccn_number) + 
      &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
                   pdq(:,:,igcm_ccn_mass) = 
@@ -1706,5 +1689,6 @@
      &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
              end where
-      END IF ! of IF (co2clouds)
+      
+      END IF                    ! of IF (co2clouds)
 
 c   9b. Aerosol particles
@@ -1750,4 +1734,8 @@
      &                tau,tauscaling)
 c Flux at the surface of co2 ice computed in co2cloud microtimestep
+           IF (co2clouds) THEN
+              zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
+           ENDIF
+
            IF (rdstorm) THEN
 c             Storm dust cannot sediment to the surface
@@ -1879,5 +1867,6 @@
           DO ig=1,ngrid
             dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
-            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
+     &                         +zdteuv(ig,l)
             pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
             pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
@@ -1894,30 +1883,17 @@
 c     (should be the last atmospherical physical process to be computed)
 c   -------------------------------------------
+
       IF (tituscap) THEN
          !!! get the actual co2 seasonal cap from Titus observations
-         CALL geticecover(ngrid, 180.*zls/pi,
+         CALL geticecover( ngrid, 180.*zls/pi,
      .                  180.*longitude/pi, 180.*latitude/pi, co2ice )
          co2ice = co2ice * 10000.
       ENDIF
       
+      
+      pdpsrf(:) = 0
 
       IF (callcond) THEN
-         zdtc(:,:) = 0.
-         zdtsurfc(:) = 0.
-         zduc(:,:) = 0.
-         zdvc(:,:) = 0.
-         zdqc(:,:,:) = 0.
-
-        if (co2clouds) then
-          CALL co2condens4micro(ngrid,nlayer,nq,ptimestep,
-     $                  capcal,zplay,zplev,tsurf,pt,
-     $                  pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
-     $                  co2ice,albedo,emis,
-     $                  zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
-     $                  fluxsurf_sw,zls,
-     $                  zdqssed_co2,zcondicea_co2microp,
-     &                  zdtcloudco2)
-        else
-          CALL co2condens(ngrid,nlayer,nq,ptimestep,
+         CALL co2condens(ngrid,nlayer,nq,ptimestep,
      $              capcal,zplay,zplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
@@ -1927,11 +1903,6 @@
      $              zdqssed_co2,zcondicea_co2microp,
      &              zdtcloudco2,zdqsc)
-           DO iq=1, nq
-           DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
-           ENDDO  ! (ig)
-           ENDDO    ! (iq)
-        end if
-        DO l=1,nlayer
+
+         DO l=1,nlayer
            DO ig=1,ngrid
              pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
@@ -1939,14 +1910,14 @@
              pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
            ENDDO
-        ENDDO
-        DO ig=1,ngrid
-           zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
-        ENDDO
-
-        IF (tracer) THEN
+         ENDDO
+         DO ig=1,ngrid
+            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+         ENDDO
+
+         IF (tracer) THEN
            DO iq=1, nq
             DO l=1,nlayer
               DO ig=1,ngrid
-                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
+                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
               ENDDO
             ENDDO
@@ -1966,4 +1937,5 @@
           ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
         ENDDO
+        
         ! update pressure levels
         DO l=1,nlayer
@@ -1974,4 +1946,5 @@
         ENDDO
         zplev(:,nlayer+1) = 0.
+        
         ! update layers altitude
         DO l=2,nlayer
@@ -1983,4 +1956,5 @@
         ENDDO
 #endif
+      
       ENDIF  ! of IF (callcond)
 
@@ -1997,4 +1971,5 @@
         ENDDO    ! (iq)
       ENDIF
+      
 c-----------------------------------------------------------------------
 c   12. Surface  and sub-surface soil temperature
@@ -2099,5 +2074,5 @@
 
       DO ig=1,ngrid
-         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
+         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig) 
       ENDDO
 
@@ -2128,5 +2103,5 @@
       ENDDO
 
-      ! Density
+      ! Density 
       DO l=1,nlayer
          DO ig=1,ngrid
@@ -2143,18 +2118,4 @@
        ENDDO
 
-      ! Total mass of co2
-      mtotco2 = 0.
-      icetotco2 = 0.
-      vaptotco2 = 0.
-      do ig=1,ngrid
-        do l=1,nlayer
-          vaptotco2 = vaptotco2 + zq(ig,l,igcm_co2) *
-     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
-          icetotco2 = icetotco2 + zq(ig,l,igcm_co2_ice) *
-     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
-        end do
-        mtotco2 = mtotco2 + co2ice(ig)
-      end do
-      mtotco2 = mtotco2 + vaptotco2 + icetotco2
 
 c    Compute surface stress : (NB: z0 is a common in surfdat.h)
@@ -2289,5 +2250,5 @@
      .          icount,' date=',ztime_fin
             
-
+            
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
@@ -2306,13 +2267,13 @@
          if (tracer) then
            ! Density-scaled opacities
-              do ig=1,ngrid
-                dsodust(ig,:) =
+              do ig=1,ngrid 
+                dsodust(ig,:) = 
      &                dsodust(ig,:)*tauscaling(ig)
-                dsords(ig,:) =
+                dsords(ig,:) = 
      &                dsords(ig,:)*tauscaling(ig)
                 dsotop(ig,:) =
      &                dsotop(ig,:)*tauscaling(ig)
               enddo
-
+         
            if(doubleq) then
               do ig=1,ngrid 
@@ -2372,4 +2333,17 @@
      &                zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
               enddo
+c D. BARDET compute integrated CO2 vapor and ice content
+              mtotco2(:)=0
+              icetotco2(:)=0
+              do ig=1,ngrid 
+               do l=1,nlayer
+                 mtotco2(ig) = mtotco2(ig) + 
+     &                      zq(ig,l,igcm_co2) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                 icetotco2(ig) = icetotco2(ig) + 
+     &                        zq(ig,l,igcm_co2_ice) * 
+     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+               enddo
+              enddo
            endif ! of if (co2clouds)
                               
@@ -2385,5 +2359,5 @@
              ENDIF !hdo
 
-             do ig=1,ngrid
+             do ig=1,ngrid 
                do l=1,nlayer
                  mtot(ig) = mtot(ig) + 
@@ -2394,9 +2368,9 @@
      &                        (zplev(ig,l) - zplev(ig,l+1)) / g
                  IF (hdo) then
-                 mtotD(ig) = mtotD(ig) +
-     &                      zq(ig,l,igcm_hdo_vap) *
+                 mtotD(ig) = mtotD(ig) + 
+     &                      zq(ig,l,igcm_hdo_vap) * 
      &                      (zplev(ig,l) - zplev(ig,l+1)) / g
-                 icetotD(ig) = icetotD(ig) +
-     &                        zq(ig,l,igcm_hdo_ice) *
+                 icetotD(ig) = icetotD(ig) + 
+     &                        zq(ig,l,igcm_hdo_ice) * 
      &                        (zplev(ig,l) - zplev(ig,l+1)) / g
                  ENDIF !hdo
@@ -2470,5 +2444,8 @@
 
            endif ! of if (water)
+
+
         endif                   ! of if (tracer)
+
 #ifndef MESOSCALE
 c        -----------------------------------------------------------------
@@ -2849,102 +2826,107 @@
 
 #else
-      !!! this is to ensure correct initialisation of mesoscale model
-      call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
-     &                 tsurf)
-      call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
-      call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
-     &                 co2ice)
-      call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
-      call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
-      call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
-      call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
-     &                 emis)
-      call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
-     &                 "K",3,tsoil)
-      call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
-     &                 "K",3,inertiedat)
+     !!! this is to ensure correct initialisation of mesoscale model
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     &                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     &                  co2ice)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     &                  emis)
+        call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
+     &                       "K",3,tsoil)
+        call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
+     &                       "K",3,inertiedat)
 #endif
+
 
 c        ----------------------------------------------------------
 c        Outputs of the CO2 cycle
 c        ----------------------------------------------------------
-
-      if (tracer.and.(igcm_co2.ne.0)) then
-!       call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
-!    &                   "kg/kg",2,zq(1,1,igcm_co2))
-        call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
-     &                   "kg/kg",3,zq(:,:,igcm_co2))
-
-        if (co2clouds) then
-          call WRITEDIAGFI(ngrid,'zdtcloudco2',
-     &                     'temperature variation of CO2 latent heat',
-     &                     'K/s',3,zdtcloudco2)
+            
+         if (tracer.and.(igcm_co2.ne.0)) then
+!          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
+!    &                     "kg/kg",2,zq(1,1,igcm_co2))
+          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
+     &                     "kg/kg",3,zq(1,1,igcm_co2))
+          if (co2clouds) then
+            CALL WRITEDIAGFI(ngrid,'mtotco2',
+     &                       'total mass of CO2 vapor',
+     &                       'kg/m2',2,mtotco2)
+            CALL WRITEDIAGFI(ngrid,'zdtcloudco2',
+     &                       'temperature variation of CO2 latent heat',
+     &                       'K/s',3,zdtcloudco2)
+
+            CALL WRITEDIAGFI(ngrid,'icetotco2',
+     &                       'total mass of CO2 ice',
+     &                       'kg/m2',2,icetotco2)
  
-          call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
-     &                     'kg/kg',3,qccnco2)
-
-          call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
-     &                     'part/kg',3,nccnco2)
-
-          call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
-     &                     3,zq(:,:,igcm_co2_ice))
-
-          call WRITEDIAGFI(ngrid,'precip_co2_ice',
-     &                     'surface deposition of co2 ice',
-     &                     'kg.m-2.s-1',2,
-     &                     zdqssed(1:ngrid,igcm_co2_ice))
-        end if ! of if (co2clouds)
-      end if ! of if (tracer.and.(igcm_co2.ne.0))
-
-      ! Output He tracer, if there is one
-      if (tracer.and.(igcm_he.ne.0)) then
-        call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
-     &                   "kg/kg",3,zq(1,1,igcm_he))
-        vmr = zq(1:ngrid,1:nlayer,igcm_he)
-     &        * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
-        call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
-     &                   'mol/mol',3,vmr)
-      end if
+            call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
+     &                       'kg/kg',3,qccnco2)
+            call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
+     &                       'part/kg',3,nccnco2)
+            call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
+     &                       3,zq(:,:,igcm_co2_ice))
+            call WRITEDIAGFI(ngrid,'precip_co2_ice',
+     &                       'surface deposition of co2 ice',
+     &                       'kg.m-2.s-1',2,
+     &                       zdqssed(1:ngrid,igcm_co2_ice))
+           endif ! of if (co2clouds)
+         endif ! of if (tracer.and.(igcm_co2.ne.0))
+        ! Output He tracer, if there is one
+        if (tracer.and.(igcm_he.ne.0)) then
+          call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
+     &                     "kg/kg",3,zq(1,1,igcm_he))
+          vmr=zq(1:ngrid,1:nlayer,igcm_he)
+     &       *mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
+          call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
+     &                       'mol/mol',3,vmr)
+        endif
 
 c        ----------------------------------------------------------
 c        Outputs of the water cycle
 c        ----------------------------------------------------------
-      if (tracer) then
-        if (water) then
+         if (tracer) then
+           if (water) then
+
 #ifdef MESOINI
             !!!! waterice = q01, voir readmeteo.F90
-          call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
-     &                     'kg/kg',3,
-     &                     zq(1:ngrid,1:nlayer,igcm_h2o_ice))
+            call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
+     &                      'kg/kg',3,
+     &                       zq(1:ngrid,1:nlayer,igcm_h2o_ice))
             !!!! watervapor = q02, voir readmeteo.F90
-          call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
-     &                     'kg/kg',3,
-     &                     zq(1:ngrid,1:nlayer,igcm_h2o_vap))
+            call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
+     &                      'kg/kg',3,
+     &                       zq(1:ngrid,1:nlayer,igcm_h2o_vap))
             !!!! surface waterice qsurf02 (voir readmeteo)
-          call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
-     &                     'kg.m-2',2,
-     &                     qsurf(1:ngrid,igcm_h2o_ice))
+            call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
+     &                      'kg.m-2',2,
+     &                       qsurf(1:ngrid,igcm_h2o_ice))
 #endif
-          call WRITEDIAGFI(ngrid,'mtot',
-     &                     'total mass of water vapor',
-     &                     'kg/m2',2,mtot)
-          call WRITEDIAGFI(ngrid,'icetot',
-     &                     'total mass of water ice',
-     &                     'kg/m2',2,icetot)
-          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
-     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
-          call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
-     &                     'mol/mol',3,vmr)
-          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
-     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
-          call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
-     &                     'mol/mol',3,vmr)
-          call WRITEDIAGFI(ngrid,'reffice',
-     &                     'Mean reff',
-     &                     'm',2,rave)
-          call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
-     &                     3,zq(:,:,igcm_h2o_ice))
-          call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
-     &                     3,zq(:,:,igcm_h2o_vap))
+            CALL WRITEDIAGFI(ngrid,'mtot',
+     &                       'total mass of water vapor',
+     &                       'kg/m2',2,mtot)
+            CALL WRITEDIAGFI(ngrid,'icetot',
+     &                       'total mass of water ice',
+     &                       'kg/m2',2,icetot)
+            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
+            CALL WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
+     &                       'mol/mol',3,vmr)
+            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
+            CALL WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
+     &                       'mol/mol',3,vmr)
+            CALL WRITEDIAGFI(ngrid,'reffice',
+     &                       'Mean reff',
+     &                       'm',2,rave)
+
+            call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
+     &             3,zq(:,:,igcm_h2o_ice))
+            call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
+     &             3,zq(:,:,igcm_h2o_vap))
 
             if (hdo) then
@@ -2972,5 +2954,5 @@
             do l=1,nlayer
                 do ig=1,ngrid
-                if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
+                if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then 
                     DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
      &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
@@ -2994,5 +2976,5 @@
             CALL WRITEDIAGFI(ngrid,'DoH_vap',
      &                       'D/H ratio in vapor',
-     &                       ' ',3,DoH_vap)
+     &                       ' ',3,DoH_vap)  
             CALL WRITEDIAGFI(ngrid,'DoH_ice',
      &                       'D/H ratio in ice',
@@ -3051,26 +3033,28 @@
               endif
 !A. Pottier
-          if (CLFvarying) then !AP14 nebulosity
-            call WRITEDIAGFI(ngrid,'totcloudfrac',
-     &                       'Total cloud fraction',
+              if (CLFvarying) then !AP14 nebulosity
+                 call WRITEDIAGFI(ngrid,'totcloudfrac',
+     &                'Total cloud fraction',
      &                       ' ',2,totcloudfrac)
-          end if !clf varying
-        end if !(water)
-
-        if (water.and..not.photochem) then
-          iq = nq
-c          write(str2(1:2),'(i2.2)') iq
-c          call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
-c    &                      'kg.m-2',2,zdqscloud(1,iq))
-c          call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
-c    &                      'kg/kg',3,zdqchim(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
-c    &                      'kg/kg',3,zdqdif(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
-c    &                      'kg/kg',3,zdqadj(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
-c    &                      'kg/kg',3,zdqc(1,1,iq))
-        end if  !(water.and..not.photochem)
-      end if !tracer
+              endif !clf varying
+
+           endif !(water)
+
+
+           if (water.and..not.photochem) then
+             iq=nq
+c            write(str2(1:2),'(i2.2)') iq
+c            call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
+c    &                       'kg.m-2',2,zdqscloud(1,iq))
+c            call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
+c    &                       'kg/kg',3,zdqchim(1,1,iq))
+c            call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
+c    &                       'kg/kg',3,zdqdif(1,1,iq))
+c            call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
+c    &                       'kg/kg',3,zdqadj(1,1,iq))
+c            call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
+c    &                       'kg/kg',3,zdqc(1,1,iq))
+           endif  !(water.and..not.photochem)
+         endif
 
 c        ----------------------------------------------------------
@@ -3078,10 +3062,10 @@
 c        ----------------------------------------------------------
 
-      call WRITEDIAGFI(ngrid,'tauref',
-     &                 'Dust ref opt depth','NU',2,tauref)
-
-      if (tracer.and.(dustbin.ne.0)) then
-
-        call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
+        call WRITEDIAGFI(ngrid,'tauref',
+     &                    'Dust ref opt depth','NU',2,tauref)
+
+         if (tracer.and.(dustbin.ne.0)) then
+
+          call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
 
 #ifndef MESOINI
@@ -3113,5 +3097,5 @@
              call WRITEDIAGFI(ngrid,'dustN','Dust number',
      &                        'part/kg',3,ndust)
-
+             
 	     select case (trim(dustiropacity))
               case ("tes")
@@ -3186,5 +3170,5 @@
              call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
      &                        'kg/kg',3,qdusttotal)
-
+             
 	     select case (trim(dustiropacity))
               case ("tes")
@@ -3211,5 +3195,5 @@
              end select
            endif ! (slpwind)
-
+	   
            if (scavenging) then
              call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
@@ -3277,5 +3261,5 @@
 c        ----------------------------------------------------------
 
-         if(callthermos) then
+         if(callthermos) then 
 
             call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
@@ -3299,13 +3283,4 @@
          endif  !(callthermos)
 
-            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
-     $           3,zdtnlte)
-            call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
-     $           3,zdteuv)
-            call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
-     $           3,zdtconduc)
-            call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
-     $           3,zdtnirco2)
-
 c        ----------------------------------------------------------
 c        ----------------------------------------------------------
@@ -3317,32 +3292,34 @@
 c        Outputs of thermals
 c        ----------------------------------------------------------
-      if (calltherm) then
+         if (calltherm) then
+
 !        call WRITEDIAGFI(ngrid,'dtke',
-!     &                   'tendance tke thermiques','m**2/s**2',
-!     &                   3,dtke_th)
+!     &              'tendance tke thermiques','m**2/s**2',
+!     &                         3,dtke_th)
 !        call WRITEDIAGFI(ngrid,'d_u_ajs',
-!     &                   'tendance u thermiques','m/s',
-!     &                   3,pdu_th*ptimestep)
+!     &              'tendance u thermiques','m/s',
+!     &                         3,pdu_th*ptimestep)
 !        call WRITEDIAGFI(ngrid,'d_v_ajs',
-!     &                   'tendance v thermiques','m/s',
-!     &                   3,pdv_th*ptimestep)
+!     &              'tendance v thermiques','m/s',
+!     &                         3,pdv_th*ptimestep)
 !        if (tracer) then
-!          if (nq .eq. 2) then
-!            call WRITEDIAGFI(ngrid,'deltaq_th',
-!     &                       'delta q thermiques','kg/kg',
-!     &                       3,ptimestep*pdq_th(:,:,2))
-!          end if
-!        end if
+!        if (nq .eq. 2) then
+!        call WRITEDIAGFI(ngrid,'deltaq_th',
+!     &              'delta q thermiques','kg/kg',
+!     &                         3,ptimestep*pdq_th(:,:,2))
+!        endif
+!        endif
 
         call WRITEDIAGFI(ngrid,'zmax_th',
-     &                   'hauteur du thermique','m',
-     &                    2,zmax_th)
+     &              'hauteur du thermique','m',
+     &                         2,zmax_th)
         call WRITEDIAGFI(ngrid,'hfmax_th',
-     &                   'maximum TH heat flux','K.m/s',
-     &                   2,hfmax_th)
+     &              'maximum TH heat flux','K.m/s',
+     &                         2,hfmax_th)
         call WRITEDIAGFI(ngrid,'wstar',
-     &                   'maximum TH vertical velocity','m/s',
-     &                   2,wstar)
-      end if
+     &              'maximum TH vertical velocity','m/s',
+     &                         2,wstar)
+
+         endif
 
 c        ----------------------------------------------------------
@@ -3382,7 +3359,8 @@
 c
 c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
-c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
-c        CALL writeg1d(ngrid,1,ps,'ps','Pa')
-c        CALL writeg1d(ngrid,nlayer,zt,'T','K')
+c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
+c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
+        
+c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
 c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
 c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
@@ -3405,7 +3383,19 @@
      &                         0,wstar)
 
-         end if ! of if (calltherm)
-
-         call WRITEDIAGFI(ngrid,'w','vertical velocity'
+         co2col(:)=0.
+         if (tracer) then
+         do l=1,nlayer
+           do ig=1,ngrid
+             co2col(ig)=co2col(ig)+
+     &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
+         enddo
+         enddo
+
+         end if
+         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
+     &                                      ,'kg/m-2',0,co2col)
+         endif ! of if (calltherm)
+
+         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
      &                              ,'m/s',1,pw)
          call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
@@ -3418,9 +3408,9 @@
          call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
          call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
-         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",
+         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
      &              "K.s-1",1,dtrad)
-         call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
+        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
      &                   'w.m-2',1,zdtsw)
-         call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
+        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
      &                   'w.m-2',1,zdtlw)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
@@ -3428,5 +3418,5 @@
 
          if (igcm_co2.ne.0) then
-          call co2sat(ngrid*nlayer,zt,zqsatco2)
+          call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
            do ig=1,ngrid 
             do l=1,nlayer
@@ -3447,5 +3437,5 @@
 c         call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
 c     &        satuco2)
-c         call WRITEDIAGFI(ngrid,"riceco2","ice radius","m"
+c         call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
 c     &        ,1,riceco2)
 ! or output in diagfi.nc (for testphys1d)
@@ -3566,7 +3556,7 @@
      &                 * (zplev(1,l) - zplev(1,l+1)) / g
                if (hdo) THEN
-                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap)
+                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap) 
      &                 * (zplev(1,l) - zplev(1,l+1)) / g
-                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice)
+                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice) 
      &                 * (zplev(1,l) - zplev(1,l+1)) / g
                ENDIF !hdo
@@ -3576,5 +3566,5 @@
                    hdotot = hdotot+mtotD+icetotD
                ENDIF ! hdo
-
+                
 
              CALL WRITEDIAGFI(ngrid,'h2otot',
@@ -3620,5 +3610,5 @@
             CALL WRITEDIAGFI(ngrid,'DoH_vap',
      &                       'D/H ratio in vapor',
-     &                       ' ',1,DoH_vap)
+     &                       ' ',1,DoH_vap)  
             CALL WRITEDIAGFI(ngrid,'DoH_ice',
      &                       'D/H ratio in ice',
@@ -3715,5 +3705,5 @@
 
         ENDIF ! hdo
-
+     
         call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
      &                    rice)
@@ -3725,28 +3715,9 @@
               endif !clfvarying
 
-        ENDIF ! of IF (water)
-
-
-      ! co2clouds
-      if (co2clouds) then
-        call WRITEDIAGFI(ngrid,'mtotco2',
-     &                   'total mass of CO2, including ice surf',
-     &                   'kg/m2',0,mtotco2)
-
-        call WRITEDIAGFI(ngrid,'vaptotco2',
-     &                   'total mass of CO2 vapor',
-     &                   'kg/m2',0,vaptotco2)
-
-        call WRITEDIAGFI(ngrid,'zdtcloudco2',
-     &                   'temperature variation of CO2 latent heat',
-     &                   'K/s',3,zdtcloudco2)
-
-        call WRITEDIAGFI(ngrid,'icetotco2',
-     &                   'total mass of CO2 ice',
-     &                   'kg/m2',0,icetotco2)
-      end if
+         ENDIF ! of IF (water)
          
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
 
          zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
@@ -3765,20 +3736,4 @@
 
       END IF       ! if(ngrid.ne.1)
-
-      co2totB = 0. ! added by C.M.
-      do ig=1,ngrid
-        do l=1,nlayer
-          co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
-     &           (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
-     &      +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
-        enddo
-        co2totB = co2totB + co2ice(ig)
-      enddo
-
-      call WRITEDIAGFI(ngrid,'co2conservation',
-     &                   'Total CO2 mass conservation in physic',
-     &                   '%',0,(co2totA-co2totB)/co2totA)
-      call test_vmr(pq(:,:,igcm_co2), pdq(:,:,igcm_co2), ptimestep,
-     &               ngrid, nlayer, 'end physiq_mod')
 
 ! XIOS outputs
Index: trunk/LMDZ.MARS/libf/phymars/tcondco2.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tcondco2.F90	(revision 2348)
+++ trunk/LMDZ.MARS/libf/phymars/tcondco2.F90	(revision 2349)
@@ -8,5 +8,5 @@
 ! Condensation temperature for co2 ice; based on 
 ! the saturation in co2sat.F JA17
-!--------------------------------------------------
+!--------------------------------------------------i
 
 integer, intent(in) :: ngrid,nlay
@@ -14,10 +14,9 @@
 double precision, intent(out), dimension(ngrid,nlay):: tcond ! CO2 condensation temperature	(atm)
 double precision:: A,B,pco2
-real :: qco2
 integer:: ig,l
 
 A=dlog(1.382d12)
 B=-3182.48
-qco2=0.
+
 DO l=1,nlay
    DO ig=1,ngrid
Index: trunk/LMDZ.MARS/libf/phymars/test_vmr_co2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/test_vmr_co2.F	(revision 2348)
+++ 	(revision )
@@ -1,22 +1,0 @@
-      subroutine test_vmr(zq, zdq, timestep, ngrid, nlay, message)
-      
-      implicit none
-      integer, intent(in)  :: ngrid, nlay
-      real, intent(in) :: zq(ngrid, nlay), zdq(ngrid, nlay), timestep
-      character(len=*), intent(in) :: message
-      integer :: i, j
-      real :: test
-      
-      do i = 1, nlay
-        do j = 1, ngrid
-          test = zq(j,i) + zdq(j,i)*timestep
-          if (test.gt.1.) then
-             write(*,*)zq(j,i), zdq(j,i), timestep
-             write(*,*)test, j, i
-             write(*,*)'Warning: mmr > 1'
-ci             call abort_physic('test_vmr', trim(message), 1)
-          end if
-        end do
-      end do
-      
-      end subroutine test_vmr
