subroutine improvedclouds(ngrid,nlay,ptimestep,
     &             pplay,pt,pdt,
     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     &             nq,tauscaling,rdust,rice,nuice,
     &             rsedcloud,rhocloud)
      implicit none
c------------------------------------------------------------------
c  This routine is used to form 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  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
c           (October 2011)
c------------------------------------------------------------------
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"
#include "comgeomfi.h"
#include "dimradmars.h"
#include "microphys.h"
c------------------------------------------------------------------
c     Inputs:

      INTEGER ngrid,nlay
      integer nq                 ! nombre de traceurs
      REAL ptimestep             ! pas de temps physique (s)
      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
      REAL pt(ngrid,nlay)        ! temperature at the middle of the
                                 !   layers (K)
      REAL pdt(ngrid,nlay)       ! tendance temperature des autres
                                 !   param.
      REAL pq(ngrid,nlay,nq)     ! traceur (kg/kg)
      REAL pdq(ngrid,nlay,nq)    ! tendance avant condensation
                                 !   (kg/kg.s-1)
      REAL tauscaling(ngridmx)     ! Convertion factor for qdust and Ndust
      REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)

c     Outputs:
      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
                                 ! (r_c in montmessin_2004)
      REAL nuice(ngrid,nlay)     ! Estimated effective variance
                                 !   of the size distribution
      REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
      REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
                                   !   H2O(kg/kg.s-1)
      REAL pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
                                   !   a la chaleur latente

c------------------------------------------------------------------
c     Local variables:

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall

      REAL*8   derf ! Error function
      !external derf

      REAL CBRT
      EXTERNAL CBRT

      INTEGER ig,l,i

      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
      REAL zt(ngridmx,nlayermx)       ! local value of temperature
      REAL zqsat(ngridmx,nlayermx)    ! saturation
      REAL lw                         !Latent heat of sublimation (J.kg-1) 
      REAL Cste
      REAL dMice           ! mass of condensated ice
      REAL sumcheck
      REAL*8 ph2o          ! Water vapor partial pressure (Pa)
      REAL*8 satu          ! Water vapor saturation ratio over ice
      REAL*8 Mo,No
      REAL*8 dN,dM,newvap
      REAL*8 Rn, Rm, dev2
      REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
      REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
      REAL*8 rate(nbin_cld)  ! nucleation rate
      REAL*8 up,dwn,Ctot,gr,seq
      REAL*8 sig      ! Water-ice/air surface tension  (N.m)
      EXTERNAL sig

c     Parameters of the size discretization
c       used by the microphysical scheme
      DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
      DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
      DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
                                           ! Minimum boundary radius (m)
      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
      DOUBLE PRECISION vrat_cld ! Volume ratio
      DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
      SAVE rb_cld
      DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m)
      DOUBLE PRECISION vol_cld(nbin_cld)   ! particle volume for each bin (m3)

      REAL sigma_ice ! Variance of the ice and CCN distributions
      SAVE sigma_ice
      
c some outputs for 1D      
      REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
      REAL dN_out(ngridmx,nlayermx)  ! mass variation for output
      REAL dM_out(ngridmx,nlayermx)  ! number variation for output
      REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
      REAL gr_out(ngridmx,nlayermx)   ! for 1d output
      REAL newvap_out(ngridmx,nlayermx)   ! for 1d output

c------------------------------------------------------------------

      IF (firstcall) THEN

c       Definition of the size grid
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~
c       rad_cld 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_cld array contains the boundary values of each rad_cld bin.
c       dr_cld is the width of each rad_cld bin.

c       Volume ratio between two adjacent bins
        vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
        vrat_cld = dexp(vrat_cld)
        write(*,*) "vrat_cld", vrat_cld

        rb_cld(1)  = rbmin_cld
        rad_cld(1) = rmin_cld
        vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.

        do i=1,nbin_cld-1
          rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
          vol_cld(i+1)  = vol_cld(i) * vrat_cld
        enddo
        
        do i=1,nbin_cld
          rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
     &      rad_cld(i)
          dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
        enddo
        rb_cld(nbin_cld+1) = rbmax_cld
        dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)

        print*, ' '
        print*,'Microphysics: size bin information:'
        print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
        print*,'-----------------------------------'
        do i=1,nbin_cld
          write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
     &      dr_cld(i)
        enddo
        write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
        print*,'-----------------------------------'

c       Contact parameter of water ice on dust ( m=cos(theta) )
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       mteta  = 0.95
        write(*,*) 'water_param contact parameter:', mteta

c       Volume of a water molecule (m3)
        vo1 = mh2o / dble(rho_ice)
c       Variance of the ice and CCN distributions
        sigma_ice = sqrt(log(1.+nuice_sed))
        
        write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
        write(*,*) 'Volume of a water molecule:', vo1

        firstcall=.false.
      END IF

c-----------------------------------------------------------------------
c     1. Initialization
c-----------------------------------------------------------------------
c     Update the needed variables

c      write(*,*) "tauscaling", tauscaling
      
      !write(*,*) "pq ccn_mass", pq(ig,:,igcm_ccn_mass)
      !write(*,*) "pdq ccn_mass", pdq(ig,:,igcm_ccn_mass)
      !write(*,*) "pq ccn_number", pq(ig,:,igcm_ccn_number)
      !write(*,*) "pdq ccn_number", pdq(ig,:,igcm_ccn_number)

!      print*, "improvedcloud debut pdq",
!     &  pdq(:,:,igcm_ccn_number)*ptimestep
!      print*, "improvedcloud debut pq", pq(:,:,igcm_ccn_number)


      do l=1,nlay
        do ig=1,ngrid
c         Temperature
          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
c         Dust mass mixing ratio
c           (converted to the true value using tauscaling)
          zq(ig,l,igcm_dust_mass) = 
     &      pq(ig,l,igcm_dust_mass) + 
     &      pdq(ig,l,igcm_dust_mass) * ptimestep
          zq(ig,l,igcm_dust_mass) = 
     &      zq(ig,l,igcm_dust_mass) * tauscaling(ig)
          zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
          zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
c         Dust particle number
c           (converted to the true value using rdust and tauscaling)
!          zq(ig,l,igcm_dust_number) = 
!     &      pq(ig,l,igcm_dust_number) + 
!     &      pdq(ig,l,igcm_dust_number) * ptimestep
          zq(ig,l,igcm_dust_number) = 
     &      (1.e0/rdust(ig,l))**3. * r3n_q * zq(ig,l,igcm_dust_mass)
          zq(ig,l,igcm_dust_number)=
     &      max(zq(ig,l,igcm_dust_number),1.E-30)
          zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
c         CCN mass mixing ratio
          zq(ig,l,igcm_ccn_mass)=
     &      pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
          zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
          zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
c          write(*,*) "pq,zq ccn_mass", pq(ig,l,igcm_ccn_mass),
c     &           zq(ig,l,igcm_ccn_mass)
c         CCN particle number
          zq(ig,l,igcm_ccn_number)=
     &      pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
          zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
          zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
c          write(*,*) "pq,zq ccn_number", pq(ig,l,igcm_ccn_number),
c     &           zq(ig,l,igcm_ccn_number)
c         Water vapor
          zq(ig,l,igcm_h2o_vap)=
     &      pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
          zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004 
          zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
c         Water ice
          zq(ig,l,igcm_h2o_ice)=
     &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 
          zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
        enddo
      enddo
      
      !print*, "improvedcloud debut pq", pq(1,:,igcm_dust_number)

c------------------------------------------------------------------
c     Cloud microphysical scheme
c------------------------------------------------------------------

      Cste = ptimestep * 4. * pi * rho_ice

      call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
            
c        write(*,*) "ccn_number avant loop phy",
c     &    zq(ig,:,igcm_ccn_number)
c        write(*,*) "ccn_mass avant loop phy",
c     &    zq(ig,:,igcm_ccn_mass)

c     Main loop over the GCM's grid
      DO l=1,nlay
      !ig = 1
        DO ig=1,ngrid

c       Get the partial pressure of water vapor and its saturation ratio
        ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
        satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
        satu_out(ig,l) = satu
        !write(*,*) "l | h2o_vap | zqsat | satu | ph2o"
        !write(*,*) l,zq(ig,l,igcm_h2o_vap), zqsat(ig,l), satu, ph2o

c       Expand the dust moments into a binned distribution
        Mo = zq(ig,l,igcm_dust_mass)
        No = zq(ig,l,igcm_dust_number)+ 1.e-30
        Rn = rdust(ig,l)
        Rm = Rn * exp( 3. * sigma_ice**2. )
        Rn = 1. / Rn
        Rm = 1. / Rm
        dev2 = 1. / ( sqrt(2.) * sigma_ice )
        do i = 1, nbin_cld
          n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 )
     &                           -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
          m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 )
     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
        enddo
        
!        sumcheck = 0
!        do i = 1, nbin_cld
!          sumcheck = sumcheck + n_aer(i)
!        enddo
!        sumcheck = abs(sumcheck/No - 1)
!        if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
!          print*, "WARNING, No sumcheck PROBLEM"
!          print*, "sumcheck, No",sumcheck, No
!          print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
!          print*, "Dust binned distribution", n_aer
!        endif
!        
!        sumcheck = 0
!        do i = 1, nbin_cld
!          sumcheck = sumcheck + M_aer(i)
!        enddo
!        sumcheck = abs(sumcheck/Mo - 1)
!        if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
!          print*, "WARNING, Mo sumcheck PROBLEM"
!          print*, "sumcheck, No",sumcheck, No
!          print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
!          print*, "Dust binned distribution", m_aer
!        endif
                
c       Get the rates of nucleation
        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)

        dN = 0.
        dM = 0.
        do i = 1, nbin_cld
          n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
          m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
          dN       = dN + n_aer(i) * rate(i) * ptimestep
          dM       = dM + m_aer(i) * rate(i) * ptimestep
        enddo
        
!        dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ), 
!     &                 zq(ig,l,igcm_dust_number) )
!     
!        dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ), 
!     &                 zq(ig,l,igcm_dust_mass) )


c        IF (zq(ig,l,igcm_ccn_number).ge.1.e-20) THEN

          Mo   = zq0(ig,l,igcm_h2o_ice) + 
     &           zq0(ig,l,igcm_ccn_mass) + 1.e-30
          No   = zq0(ig,l,igcm_ccn_number)
          !write(*,*) "l,cloud particles,cloud mass",l, No, Mo
          rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
     &                     +zq0(ig,l,igcm_ccn_mass) / Mo * rho_dust
          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
          rice(ig,l) =
     &      ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
          if (Mo.lt.1.e-20) rice(ig,l) = 1.e-8
          seq  = exp( 2.*sig(zt(ig,l))*mh2o / 
     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )

          call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
     &                    ph2o,ph2o/satu,seq,rice(ig,l),gr)

          up  = Cste * gr * rice(ig,l) * No * seq + 
     &            zq(ig,l,igcm_h2o_vap)
          dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
  
          Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
          
          newvap = min(up/dwn,Ctot)
  
          gr = gr * ( newvap/zqsat(ig,l) - seq )

          
          dMice = min( max(Cste * No * rice(ig,l) * gr,
     &                     -zq(ig,l,igcm_h2o_ice) ), 
     &                 zq(ig,l,igcm_h2o_vap) )
     
c----------- TESTS 1D output ---------
          if (ngrid.eq.1) then
             Mcon_out(ig,l) = dMice
             newvap_out(ig,l) = newvap
             gr_out(ig,l) = gr
             dN_out(ig,l) = dN
             dM_out(ig,l) = dM
          endif
c-------------------------------------
          
c       Water ice
          zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) + 
     &        dMice
     
c       Water vapor
          zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - 
     &        dMice

  
c         If all the ice particles sublimate, all the condensation
c           nuclei are released:
          if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
c           Water ice particles
            zq(ig,l,igcm_h2o_ice) = 0.
c           Dust particles
            zq(ig,l,igcm_dust_mass ) = zq0(ig,l,igcm_dust_mass) + 
     &        zq0(ig,l,igcm_ccn_mass)
            zq(ig,l,igcm_dust_number ) = zq0(ig,l,igcm_dust_number) + 
     &        zq0(ig,l,igcm_ccn_number)
c           CCNs
            zq(ig,l,igcm_ccn_mass) = 0.
            zq(ig,l,igcm_ccn_number) = 0.
c           for coherence 
            dM = 0
            dN = 0
          endif
c        ELSE
cc         Initialize rhocloud and rice to avoid divisions by 0
c          rhocloud(ig,l) = 1.e-10
c          rice(ig,l) = 1.e-8
c          dM = 0
c          dN = 0
c        ENDIF ! of if (zq(ig,l,igcm_ccn_number).ge.1.e-20)

c       Dust particles
        zq(ig,l,igcm_dust_mass ) = zq(ig,l,igcm_dust_mass ) - dM
        zq(ig,l,igcm_dust_number ) = zq(ig,l,igcm_dust_number ) - dN
c       CCNs
        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN

        ENDDO
      ENDDO


!      print*, "improvedclouds zq0 abs.", zq0(:,:,igcm_ccn_number)
!      print*, "improvedclouds zq abs.", zq(:,:,igcm_ccn_number)
c------------------------------------------------------------------
c     Convert the initial values back into relative values
c       (has to be done before updating rdust!)
c------------------------------------------------------------------

      do l=1, nlay
        do ig=1,ngrid
          zq0(ig,l,igcm_dust_mass) =
     &      zq0(ig,l,igcm_dust_mass) / tauscaling(ig)
          zq0(ig,l,igcm_dust_number) =
     &      (1.e0/rdust(ig,l))**3. * r3n_q * 
     &      zq0(ig,l,igcm_dust_mass)
        enddo ! of do ig=1,ngrid
      enddo ! of do l=1,nlay

c------------------------------------------------------------------
c     Update the dust radius
c------------------------------------------------------------------

      DO l=1,nlay
        DO ig=1,ngrid
          rdust(ig,l)=
     &      CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
     &      max(zq(ig,l,igcm_dust_number),0.01))
          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
        ENDDO
      ENDDO

c------------------------------------------------------------------
c     Convert zq back into relative values (only applies to dust)
c------------------------------------------------------------------

      do l=1, nlay
        do ig=1,ngrid
c         Dust mass mixing ratio
c           (converted back into relative value using tauscaling)
          zq(ig,l,igcm_dust_mass) =
     &      zq(ig,l,igcm_dust_mass) / tauscaling(ig)
          zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
c         Dust particle number
c           (converted back into relative value)
          zq(ig,l,igcm_dust_number) =
     &      (1.e0/rdust(ig,l))**3. * r3n_q * zq(ig,l,igcm_dust_mass)
          zq(ig,l,igcm_dust_number)=
     &      max(zq(ig,l,igcm_dust_number),1.E-30)
        enddo ! of do ig=1,ngrid
      enddo ! of do l=1,nlay
      
!      print*, "improvedclouds zq0 rel.", zq0(1,:,igcm_ccn_number)
!      print*, "improvedclouds zq rel.", zq(1,:,igcm_ccn_number)
      
c------------------------------------------------------------------
c     Compute the sedimentation radius
c------------------------------------------------------------------

      do l=1, nlay
        do ig=1,ngrid
          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
     &                         rdust(ig,l) )
          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
        enddo ! of do ig=1,ngrid
      enddo ! of do l=1,nlay

c------------------------------------------------------------------
c     Force positive values
c------------------------------------------------------------------

!      do l=1, nlay
!        do ig=1,ngrid
!          zq(ig,l,igcm_ccn_mass)=
!     &       max(zq(ig,l,igcm_ccn_mass),1e-30)
!          zq(ig,l,igcm_ccn_number)=
!     &       max(zq(ig,l,igcm_ccn_number),1e-30)
!          zq(ig,l,igcm_h2o_vap)=
!     &       max(zq(ig,l,igcm_h2o_vap),1e-30)
!          zq(ig,l,igcm_h2o_ice)=
!     &       max(zq(ig,l,igcm_h2o_ice),1e-30)
!        enddo ! of do ig=1,ngrid
!      enddo ! of do l=1,nlay


c------------------------------------------------------------------
c     Compute the final tendencies
c------------------------------------------------------------------

c     Initialize the tendencies
      pdqscloud(1:ngrid,1:nq)=0
      pdqcloud(1:ngrid,1:nlay,1:nq)=0
      pdtcloud(1:ngrid,1:nlay)=0

c     Update the tendencies
      do l=1, nlay
        do ig=1,ngrid
          pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
     &                            -zq0(ig,l,igcm_dust_mass))/ptimestep
          pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
     &                            -pq(ig,l,igcm_dust_number))/ptimestep
     &                            - pdq(ig,l,igcm_dust_number) !!! AJOUT TN
          pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
     &                            -zq0(ig,l,igcm_ccn_mass))/ptimestep
          pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
     &                           -zq0(ig,l,igcm_ccn_number))/ptimestep
          pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
     &                            -zq0(ig,l,igcm_h2o_vap))/ptimestep
          pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
     &                            -zq0(ig,l,igcm_h2o_ice))/ptimestep
          lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
          pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
        end do
      end do
!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
!     &                    satu_out)

!      print*, "improvedclouds pdq*dt", 
!     &         pdqcloud(:,:,igcm_ccn_number)*ptimestep


c------------------------------------------------------------------
c     TEST_JBM
      IF (ngrid.eq.1) THEN
c         call WRITEDIAGFI(ngrid,"tausca","tauscaling","",0,
c     &                    tauscaling)
         call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
     &                    newvap_out)
         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
     &                    gr_out)
         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
     &                    dM_out)
         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
     &                    dN_out)
         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
     &                    Mcon_out)
         call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
     &                    zqsat)
         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
     &                    satu_out)
         call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
     &                    rice)
         call WRITEDIAGFI(ngrid,"rdust","rdust","m",1,
     &                    rdust)
         call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
     &                    rsedcloud)
         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
     &                    rhocloud)
      ENDIF
c------------------------------------------------------------------
      return
      end