SUBROUTINE callsedim_pluto(ngrid,nlay, ptimestep, &
                     pplev,zlev,pt,pdt,rice_ch4,rice_co, &
                     pq, pdqfi, pdqsed,pdqs_sed,nq,pphi)

      use radinc_h, only : naerkind
      use tracer_h, only: igcm_ch4_ice,igcm_co_ice,radius,rho_q
      use comcstfi_mod, only:  g
      IMPLICIT NONE
!==================================================================
!
!     Purpose
!     -------
!     Calculates sedimentation of aerosols depending on their
!     density and radius.
!
!     Authors
!     -------
!     F. Forget (1999)
!     Tracer generalisation by E. Millour (2009)
!
!==================================================================

!-----------------------------------------------------------------------
!   declarations:
!   -------------

!
!   arguments:
!   ----------

      INTEGER ngrid              ! number of horizontal grid points
      INTEGER nlay               ! number of atmospheric layers
      REAL ptimestep             ! physics time step (s)
      REAL pplev(ngrid,nlay+1)   ! pressure at inter-layers (Pa)
      REAL pt(ngrid,nlay)        ! temperature at mid-layer (K)
      REAL pdt(ngrid,nlay)       ! tendency on temperature
      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
      REAL pphi(ngrid,nlay)      ! geopotential


!    Traceurs :
      integer nq             ! number of tracers
      real pq(ngrid,nlay,nq)  ! tracers (kg/kg)
      real pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
      real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
      real pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)

!   local:
!   ------

      INTEGER l,ig, iq

      real zqi(ngrid,nlay) ! to locally store tracers
      real zt(ngrid,nlay) ! to locally store temperature (K)
      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
      real epaisseur (ngrid,nlay) ! Layer thickness (m)
      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
      real rfall_ch4(ngrid,nlay)
      real rfall_co(ngrid,nlay)
      real rice_ch4(ngrid,nlay)
      real rice_co(ngrid,nlay)

      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./

!    ** un petit test de coherence
!       --------------------------

      IF (firstcall) THEN
         IF(ngrid.NE.ngrid) THEN
            PRINT*,'STOP dans callsedim'
            PRINT*,'probleme de dimensions :'
            PRINT*,'ngrid  =',ngrid
            PRINT*,'ngrid  =',ngrid
            STOP
         ENDIF

        firstcall=.false.
      ENDIF ! of IF (firstcall)

!=======================================================================
!     Preliminary calculation of the layer characteristics
!     (mass (kg.m-2), thickness (m), etc.

      do  l=1,nlay
        do ig=1, ngrid
          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
        end do
      end do

      do iq=1,nq
        if(radius(iq).gt.1.e-12) then   ! no sedimentation for gases (defined by radius=0)
!     Radius values are defined in initracer
!     The value of q is updated after the other parameterisations

          do l=1,nlay
            do ig=1,ngrid
              ! store locally updated tracers
              zqi(ig,l)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep

!                cf sur Mars:
!                On affecte un rayon moyen aux poussieres a chaque altitude du type :
!                r(z)=r0*exp(-z/H) avec r0=0.8 micron et H=18 km.
!                rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )
!                Pluton : choix de rfall a la place de radius
              if (iq.eq.igcm_ch4_ice) then
                 ! TB: rfall_ch4(ig,l)=max( rice_ch4(ig,l)*1.5,2.e-7)
                 rfall_ch4(ig,l)=max( rice_ch4(ig,l),2.e-7)
                 rfall_ch4(ig,l)=min(rfall_ch4(ig,l),1.e-4)
              endif
              if (iq.eq.igcm_co_ice) then
                 rfall_co(ig,l)=max( rice_co(ig,l),2.e-7)
                 rfall_co(ig,l)=min(rfall_co(ig,l),1.e-4)
              endif
            enddo
          enddo ! of do l=1,nlay

!=======================================================================
!     Calculate the transport due to sedimentation for each tracer

          if (iq.eq.igcm_ch4_ice) then
          !if (iceparty.and.(iq.eq.igcm_ch4_ice)) then
            call newsedim_pluto(ngrid,nlay,ngrid*nlay,ptimestep, &
           pplev,masse,epaisseur,zt,rfall_ch4,rho_q(iq),zqi,wq,pphi)
          else if (iq.eq.igcm_co_ice) then
            call newsedim_pluto(ngrid,nlay,ngrid*nlay,ptimestep, &
           pplev,masse,epaisseur,zt,rfall_co,rho_q(iq),zqi,wq,pphi)
          else if ((radius(iq).gt.0.)) then   ! haze tracers
            call newsedim_pluto(ngrid,nlay,1,ptimestep, &
           pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),zqi,wq,pphi)
          endif

!=======================================================================
!     Calculate the tendencies

          do ig=1,ngrid
!     Ehouarn: with new way of tracking tracers by name, this is simply
            pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
          end do

          DO l = 1, nlay
            DO ig=1,ngrid
              pdqsed(ig,l,iq)=(zqi(ig,l)- &
             (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
            ENDDO
          ENDDO

        endif ! of if(radius(iq).gt.1.e-12)
      enddo ! of do iq=1,nq

      RETURN
      END