SUBROUTINE callsedim(ngrid,nlay, ptimestep,
     &                pplev,zlev, zlay, pt, rdust, rice,
     &                rsedcloud,rhocloud,
     &                pq, pdqfi, pdqsed,pdqs_sed,nq, 
     &                tau, tauscaling)
      IMPLICIT NONE

c=======================================================================
c      Sedimentation of the  Martian aerosols
c      depending on their density and radius
c
c      F.Forget 1999
c
c      Modified by J.-B. Madeleine 2010: Now includes the doubleq
c        technique in order to have only one call to callsedim in
c        physiq.F.
c
c=======================================================================

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

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

c
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    Aerosol radius provided by the water ice microphysical scheme:
      REAL rdust(ngrid,nlay)     ! Dust geometric mean radius (m)
      REAL rice(ngrid,nlay)      ! Ice geometric mean radius (m)

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

      REAL CBRT
      EXTERNAL CBRT

      INTEGER l,ig, iq

      real zqi(ngridmx,nlayermx,nqmx) ! 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)
      real r0(ngridmx,nlayermx) ! geometric mean radius used for
                                !   sedimentation (m)
      real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
                                    !   dust (m)
      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
                                    !   CCNs (m)
c     Sedimentation radius of water ice
      real rsedcloud(ngridmx,nlayermx)
c     Cloud density (kg.m-3)
      real rhocloud(ngridmx,nlayermx)
      
c     for ice radius computation
      REAL ccn_factor
      REAL Mo,No
      REAL tau(ngrid,nlay), tauscaling(ngrid)
      REAL zlay(ngrid,nlay)   ! altitude at the middle of the layers
      REAl ccntyp



c     Discrete size distributions (doubleq)
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c       1) Parameters used to represent the changes in fall
c          velocity as a function of particle size;
      integer nr,ir
      parameter (nr=12) !(nr=7) ! number of bins
      real rd(nr),qr(ngridmx,nlayermx,nr)
      real rdi(nr+1)    ! extreme and intermediate radii
      real Sq(ngridmx,nlayermx)
      real rdmin,rdmax,rdimin,rdimax
      data rdmin/1.e-8/ !/1.e-7/
      data rdmax/30.e-6/
      data rdimin/1.e-10/
      data rdimax/1e-4/
      save rd, rdi

c       2) Second size distribution for the log-normal integration
c          (the mass mixing ratio is computed for each radius)

      integer ninter, iint
      parameter (ninter=4)  ! nombre de point entre chaque rayon rdi
      real rr(ninter,nr)
      save rr
      integer radpower
      real sigma0

c       3) Other local variables used in doubleq

      INTEGER idust_mass   ! index of tracer containing dust mass
                           !   mix. ratio
      INTEGER idust_number ! index of tracer containing dust number
                           !   mix. ratio
      INTEGER iccn_mass    ! index of tracer containing CCN mass
                           !   mix. ratio
      INTEGER iccn_number  ! index of tracer containing CCN number
                           !   mix. ratio
      SAVE idust_mass,idust_number
      SAVE iccn_mass,iccn_number

c     Firstcall:

      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

c       Doubleq: initialization
        IF (doubleq) THEN
         do ir=1,nr
             rd(ir)= rdmin*(rdmax/rdmin)**(float(ir-1)/float(nr-1))
         end do
         rdi(1)=rdimin
         do ir=2,nr
           rdi(ir)= sqrt(rd(ir-1)*rd(ir))
         end do
         rdi(nr+1)=rdimax

         do ir=1,nr
           do iint=1,ninter
             rr(iint,ir)=
     &        rdi(ir)*
     &        (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
c             write(*,*) rr(iint,ir)
           end do
         end do

      ! identify tracers corresponding to mass mixing ratio and
      ! number mixing ratio
        idust_mass=0      ! dummy initialization
        idust_number=0    ! dummy initialization

        do iq=1,nq
          if (noms(iq).eq."dust_mass") then
            idust_mass=iq
          endif
          if (noms(iq).eq."dust_number") then
            idust_number=iq
          endif
        enddo

        ! check that we did find the tracers
        if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
          write(*,*) 'callsedim: error! could not identify'
          write(*,*) ' tracers for dust mass and number mixing'
          write(*,*) ' ratio and doubleq is activated!'
          stop
        endif
        ENDIF !of if (doubleq)

        IF (scavenging) THEN
          iccn_mass=0
          iccn_number=0
          do iq=1,nq
            if (noms(iq).eq."ccn_mass") then
              iccn_mass=iq
            endif
            if (noms(iq).eq."ccn_number") then
              iccn_number=iq
            endif
          enddo
          ! check that we did find the tracers
          if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
            write(*,*) 'callsedim: error! could not identify'
            write(*,*) ' tracers for ccn mass and number mixing'
            write(*,*) ' ratio and scavenging is activated!'
            stop
          endif
        ELSE
          write(*,*) "water_param CCN reduc. fac. ", ccn_factor
          write(*,*) "Careful: only used when microphys=F, otherwise"
          write(*,*) "  the contact parameter is used instead;"
        ENDIF !of if (scavenging)

        IF (water) THEN
          write(*,*) "water_param nueff Sedimentation:", nuice_sed
          IF (activice) THEN
            write(*,*) "water_param nueff Radiative:", nuice_ref
          ENDIF
        ENDIF
      
        firstcall=.false.
      ENDIF ! of IF (firstcall)

c-----------------------------------------------------------------------
c    1. Initialization
c    -----------------

      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
c     Updating the mass mixing ratio with the tendencies coming
c       from other parameterizations:
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      do iq=1,nq
        do l=1,nlay
          do ig=1,ngrid
            zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
          enddo
        enddo 
      enddo

c    Computing the different layer properties
c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c    Mass (kg.m-2), thickness(m), crossing time (s)  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

c =================================================================
c     Compute the geometric mean radius used for sedimentation

      if (doubleq) then
        do l=1,nlay
          do ig=1, ngrid
            r0dust(ig,l) =
     &        CBRT(r3n_q*zqi(ig,l,idust_mass)/
     &        max(zqi(ig,l,idust_number),0.01))
            r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
          end do
        end do
      endif
      if (scavenging) then
        do l=1,nlay
          do ig=1, ngrid
            r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
          end do
        end do
      endif

c =================================================================
      do iq=1,nq
        if(radius(iq).gt.1.e-9) then   ! no sedim for gaz

c -----------------------------------------------------------------
c         DOUBLEQ CASE
c -----------------------------------------------------------------
          if ((doubleq.and.
     &        ((iq.eq.idust_mass).or.
     &         (iq.eq.idust_number))).or.
     &        (scavenging.and.
     &        ((iq.eq.iccn_mass).or.
     &        (iq.eq.iccn_number)))) then

c           Computing size distribution:
c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
              do  l=1,nlay
                do ig=1, ngrid
                  r0(ig,l)=r0dust(ig,l)
                end do
              end do
              sigma0 = varian
            else 
              do  l=1,nlay
                do ig=1, ngrid
                  r0(ig,l)=r0ccn(ig,l)
                end do
              end do
              sigma0 = sqrt(log(1.+nuice_sed))
            endif

c        Computing mass mixing ratio for each particle size
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
            radpower = 2
          ELSE
            radpower = -1
          ENDIF
          Sq(1:ngrid,1:nlay) = 0.
          do ir=1,nr
            do l=1,nlay
              do ig=1,ngrid
c                ****************
c                Size distribution integration
c                (Trapezoid Integration Method)
                 qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))*
     &             (rr(1,ir)**radpower)*
     &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*sigma0**2))
                 do iint=2,ninter-1
                   qr(ig,l,ir)=qr(ig,l,ir) +
     &             0.5*(rr(iint+1,ir)-rr(iint-1,ir))*
     &             (rr(iint,ir)**radpower)*
     &             exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
     &             (2*sigma0**2))
                 end do
                 qr(ig,l,ir)=qr(ig,l,ir) +
     &             0.5*(rr(ninter,ir)-rr(ninter-1,ir))*
     &             (rr(ninter,ir)**radpower)*
     &             exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
     &             (2*sigma0**2))

c                **************** old method (not recommended!)
c                qr(ig,l,ir)=(rd(ir)**(5-3*iq))*
c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*sigma0**2) )
c                ******************************

                 Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
              enddo
            enddo
          enddo

          do ir=1,nr
            do l=1,nlay
              do ig=1,ngrid
                 qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
              enddo
            enddo
          enddo

c         Computing sedimentation for each tracer
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

          zqi(1:ngrid,1:nlay,iq) = 0.
          pdqs_sed(1:ngrid,iq) = 0.

          do ir=1,nr
             IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
               call newsedim(ngrid,nlay,1,1,ptimestep,
     &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
     &         wq,0.5)
             ELSE
               call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
     &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
     &         wq,1.0)
             ENDIF

c            Tendencies
c            ~~~~~~~~~~
             do ig=1,ngrid
               pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
     &                                + wq(ig,1)/ptimestep
             end do
             DO l = 1, nlay
               DO ig=1,ngrid
                 zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
               ENDDO
             ENDDO
          enddo ! of do ir=1,nr
c -----------------------------------------------------------------
c         WATER CYCLE CASE
c -----------------------------------------------------------------
          else if (water.and.(iq.eq.igcm_h2o_ice)) then
            if (microphys) then
              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
     &        zqi(1,1,iq),wq,1.0)
            else
              call newsedim(ngrid,nlay,ngrid*nlay,1,
     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
     &        zqi(1,1,iq),wq,1.0)
            endif ! of if (microphys)
c           Tendencies
c           ~~~~~~~~~~
            do ig=1,ngrid 
              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
            end do
c -----------------------------------------------------------------
c         GENERAL CASE
c -----------------------------------------------------------------
          else
            call newsedim(ngrid,nlay,1,1,ptimestep,
     &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
     &      zqi(1,1,iq),wq,1.0)
c           Tendencies
c           ~~~~~~~~~~
            do ig=1,ngrid 
              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
            end do
          endif ! of if doubleq and if water
c -----------------------------------------------------------------

c         Compute the final tendency:
c         ---------------------------
          DO l = 1, nlay
            DO ig=1,ngrid
              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
     $        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
            ENDDO
          ENDDO

        endif ! of if(radius(iq).gt.1.e-9)
c =================================================================
      enddo ! of do iq=1,nq
 
c     Update the dust particle size "rdust"
c     -------------------------------------
      DO l = 1, nlay
        DO ig=1,ngrid
          rdust(ig,l)=
     &      CBRT(r3n_q*zqi(ig,l,idust_mass)/
     &      max(zqi(ig,l,idust_number),0.01))
          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
        ENDDO
      ENDDO
      
c     Update the ice particle size "rice"
c     -------------------------------------
      IF(scavenging) THEN 
        DO l = 1, nlay
          DO ig=1,ngrid
            Mo   = zqi(ig,l,igcm_h2o_ice) + 
     &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
            No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
            rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
     &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
     &                        / 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.)
           if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
c           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
c           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
           
          ENDDO
        ENDDO
      ELSE
        DO l = 1, nlay
          DO ig=1,ngrid
            ccntyp = 
     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
            ccntyp = ccntyp /ccn_factor
            rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))      
          ENDDO
        ENDDO
      ENDIF

      RETURN
      END