SUBROUTINE callsedim(ngrid,nlay, ptimestep,
     $                pplev,zlev, pt,
     &                pq, pdqfi, pdqsed,pdqs_sed,nq,rfall)
      IMPLICIT NONE

!==================================================================
!     
!     Purpose
!     -------
!     Calculates sedimentation of aerosols depending on their
!     density and radius.
!     
!     Authors
!     -------
!     F. Forget (1999)
!     Tracer generalisation by E. Millour (2009)
!     
!==================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "tracer.h"
#include "callkeys.h"

#include "fisice.h"

c   arguments:
c   ----------

      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 zlev(ngrid,nlay+1)    ! altitude at layer boundaries


c    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)
      
c   local:
c   ------

      INTEGER l,ig, iq

      real zqi(ngridmx,nlayermx) ! to locally store tracers
      real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
      real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
      real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
c      real dens(ngridmx,nlayermx) ! Mean density of the ice part. accounting for dust core
      real rfall(ngridmx,nlayermx)


      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./

c    ** un petit test de coherence
c       --------------------------

      IF (firstcall) THEN
         IF(ngrid.NE.ngridmx) THEN
            PRINT*,'STOP dans callsedim'
            PRINT*,'probleme de dimensions :'
            PRINT*,'ngrid  =',ngrid
            PRINT*,'ngridmx  =',ngridmx
            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)
        end do
      end do

      iq=igcm_h2o_ice

!     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
            enddo
          enddo ! of do l=1,nlay



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

          call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
     &         pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq)

!=======================================================================
!     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

      return
      end