!----- This subroutine calculates the sedimentation flux of Tracers
!
SUBROUTINE sediment_mod(t_seri,pplay,zrho,paprs,time_step,RHcl, &
        id_coss,id_codu,id_scdu, &
        ok_chimeredust, &
        sed_ss,sed_dust,sed_dustsco, &
        sed_ss3D,sed_dust3D,sed_dustsco3D,tr_seri)
  !nhl     .                                       xlon,xlat,
  !
USE chem_mod_h
  USE yoecumf_mod_h
     USE dimphy
  USE infotrac_phy, ONLY: nbtr
USE yomcst_mod_h
IMPLICIT NONE
  !

    ! INCLUDE "dimphy.h"

  !
   REAL :: RHcl(klon,klev)     ! humidite relative ciel clair
   REAL :: tr_seri(klon, klev,nbtr) !conc of tracers
   REAL :: sed_ss(klon) !sedimentation flux of Sea Salt (g/m2/s)
   REAL :: sed_dust(klon) !sedimentation flux of dust (g/m2/s)
   REAL :: sed_dustsco(klon) !sedimentation flux of scoarse  dust (g/m2/s)
   REAL :: sed_ss3D(klon,klev) !sedimentation flux of Sea Salt (g/m2/s)
   REAL :: sed_dust3D(klon,klev) !sedimentation flux of dust (g/m2/s)
   REAL :: sed_dustsco3D(klon,klev) !sedimentation flux of scoarse  dust (g/m2/s)
   REAL :: t_seri(klon, klev)   !Temperature at mid points of Z (K)
   REAL :: v_dep_ss(klon,klev)  ! sed. velocity for SS m/s
   REAL :: v_dep_dust(klon,klev)  ! sed. velocity for dust m/s
   REAL :: v_dep_dustsco(klon,klev)  ! sed. velocity for dust m/s
   REAL :: pplay(klon, klev)    !pressure at mid points of Z (Pa)
   REAL :: zrho(klon, klev)     !Density of air at mid points of Z (kg/m3)
   REAL :: paprs(klon, klev+1)    !pressure at interface of layers Z (Pa)
   REAL :: time_step            !time step (sec)
   LOGICAL :: ok_chimeredust
   REAL :: xlat(klon)       ! latitudes pour chaque point
   REAL :: xlon(klon)       ! longitudes pour chaque point
   INTEGER :: id_coss,id_codu,id_scdu
  !
  !------local variables
  !
   INTEGER :: i, k, nbre_RH
   PARAMETER(nbre_RH=12)
  !
   REAL :: lambda, ss_g
   REAL :: mmd_ss      !mass median diameter of SS (um)
   REAL :: mmd_dust          !mass median diameter of dust (um)
   REAL :: mmd_dustsco          !mass median diameter of scoarse dust (um)
   REAL :: rho_ss(nbre_RH),rho_ss1 !density of sea salt (kg/m3)
   REAL :: rho_dust          !density of dust(kg/m3)
   REAL :: v_stokes, CC, v_sed, ss_growth_f(nbre_RH)
   REAL :: sed_flux(klon,klev)  ! sedimentation flux g/m2/s
   REAL :: air_visco(klon,klev)
   REAL :: zdz(klon,klev)       ! layers height (m)
   REAL :: temp                 ! temperature in degree Celius
  !
   INTEGER :: RH_num
   REAL :: RH_MAX, DELTA, rh, RH_tab(nbre_RH)
   PARAMETER (RH_MAX=95.)
  !
   DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
  !
  !
   DATA rho_ss/2160. ,2160. ,2160.,  2160,  1451.6, 1367.9, &
         1302.9,1243.2,1182.7, 1149.5,1111.6, 1063.1/
  !
   DATA ss_growth_f/0.503, 0.503, 0.503, 0.503, 0.724, 0.782, &
         0.838, 0.905, 1.000, 1.072, 1.188, 1.447/
  !
  !
   mmd_ss=12.7   !dia -um at 80% for bin 0.5-20 um but 90% of real mmd
            ! obsolete      mmd_dust=2.8  !micrometer for bin 0.5-20 and 0.5-10 um
  ! 4tracer SPLA:       mmd_dust=11.0  !micrometer for bin 0.5-20 and 0.5-10 um
  !3days       mmd_dust=3.333464  !micrometer for bin 0.5-20 and 0.5-10 um
  !3days       mmd_dustsco=12.91315  !micrometer for bin 0.5-20 and 0.5-10 um
  !JE20140911       mmd_dust=3.002283  !micrometer for bin 0.5-20 and 0.5-10 um
  !JE20140911       mmd_dustsco=13.09771  !micrometer for bin 0.5-20 and 0.5-10 um
  !JE20140911        mmd_dust=5.156346  !micrometer for bin 0.5-20 and 0.5-10 um
  !JE20140911        mmd_dustsco=15.56554  !micrometer for bin 0.5-20 and 0.5-10 um
    IF (ok_chimeredust) THEN
  !JE20150212<< : changes in ustar in dustmod changes emission distribution
     ! mmd_dust=3.761212  !micrometer for bin 0.5-3 and 0.5-10 um
     ! mmd_dustsco=15.06167  !micrometer for bin 3-20 and 0.5-10 um
  !JE20150212>>
  !JE20150618: Change in div3 of dustmod changes distribution. now is div3=6
  !div=3        mmd_dust=3.983763
  !div=3        mmd_dustsco=15.10854
    mmd_dust=3.898047
    mmd_dustsco=15.06167
    ELSE
    mmd_dust=11.0  !micrometer for bin 0.5-20 and 0.5-10 um
    mmd_dustsco=100. ! absurd value, bin not used in this scheme
    ENDIF


   rho_dust=2600. !kg/m3
  !
  !--------- Air viscosity (poise=0.1 kg/m-sec)-----------
  !
   DO k=1, klev
   DO i=1, klon
  !
   zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
  !
   temp=t_seri(i,k)-RTT
  !
   IF (temp.LT.0.) THEN
     air_visco(i,k)=(1.718+0.0049*temp-1.2e-5*temp*temp)*1.e-4
   ELSE
     air_visco(i,k)=(1.718+0.0049*temp)*1.e-4
   ENDIF
  !
   ENDDO
   ENDDO
  !
  !--------- for Sea Salt -------------------
  !
  !
  !
   IF(id_coss>0) THEN
   DO k=1, klev
   DO i=1,klon
  !
  !---cal. correction factor hygroscopic growth of aerosols
  !
    rh=MIN(RHcl(i,k)*100.,RH_MAX)
    RH_num = INT( rh/10. + 1.)
    IF (rh.gt.85.) RH_num=10
    IF (rh.gt.90.) RH_num=11
    DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
  !
    ss_g=ss_growth_f(rh_num) + &
          DELTA*(ss_growth_f(RH_num+1)-ss_growth_f(RH_num))

    rho_ss1=rho_ss(rh_num) + &
          DELTA*(rho_ss(RH_num+1)-rho_ss(RH_num))
  !
    v_stokes=RG*(rho_ss1-zrho(i,k))* & !m/sec
          (mmd_ss*ss_g)*(mmd_ss*ss_g)* &
          1.e-12/(18.0*air_visco(i,k)/10.)
  !
   lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
  !
   CC=1.0+1.257*lambda/(mmd_ss*ss_g)/1.e6  ! C-correction factor
  !
   v_sed=v_stokes*CC                       ! m/sec !orig
  !
  !---------check for v_sed*dt<zdz
  !
   IF (v_sed*time_step.GT.zdz(i,k)) THEN
     v_sed=zdz(i,k)/time_step
   ENDIF
  !
   v_dep_ss(i,k)= v_sed
   sed_flux(i,k)= tr_seri(i,k,id_coss)*v_sed !g/cm3*m/sec
   ! !sed_ss3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
  ! ! conc_sed_ss3D(i,k)=sed_flux(i,k)*1.e6      !g/m3*sec !!!!!!!
  !
   ENDDO          !klon
   ENDDO          !klev
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   sed_ss3D(:,:)=0.0  ! initialisation

   DO k=1, klev
   DO i=1, klon
   sed_ss3D(i,k)=sed_ss3D(i,k)- &
         sed_flux(i,k)/zdz(i,k) !!!!!!!!!!!!!!!!!!!!!!
   ENDDO          !klon
   ENDDO          !klev
  !
   DO k=1, klev-1
   DO i=1, klon
    sed_ss3D(i,k)=sed_ss3D(i,k)+ &
          sed_flux(i,k+1)/zdz(i,k) !!!!!!!!

   ENDDO          !klon
   ENDDO          !klev

  DO k = 1, klev
  DO i = 1, klon
      tr_seri(i,k,id_coss)=tr_seri(i,k,id_coss)+ &
            sed_ss3D(i,k)*time_step
  ENDDO
  ENDDO

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
   DO i=1, klon
     sed_ss(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
   ENDDO          !klon
   ELSE
    DO i=1, klon
      sed_ss(i)=0.
    ENDDO
   ENDIF
  !
  !

  !--------- For dust ------------------
  !
  !
   IF(id_codu>0) THEN
   DO k=1, klev
   DO i=1,klon
  !
    v_stokes=RG*(rho_dust-zrho(i,k))* & !m/sec
          mmd_dust*mmd_dust* &
          1.e-12/(18.0*air_visco(i,k)/10.)
  !
   lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
   CC=1.0+1.257*lambda/(mmd_dust)/1.e6        !dimensionless
   v_sed=v_stokes*CC                       !m/sec
  !
  !---------check for v_sed*dt<zdz
  !
   IF (v_sed*time_step.GT.zdz(i,k)) THEN
     v_sed=zdz(i,k)/time_step
   ENDIF

  !
   v_dep_dust(i,k)= v_sed
   sed_flux(i,k)  = tr_seri(i,k,id_codu)*v_sed !g/cm3.m/sec
   ! !sed_dust3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
  !
   ENDDO          !klon
   ENDDO          !klev

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   sed_dust3D(:,:)=0.0  ! initialisation

   DO k=1, klev
   DO i=1, klon
   sed_dust3D(i,k)=sed_dust3D(i,k)- &
         sed_flux(i,k)/zdz(i,k)
   ENDDO          !klon
   ENDDO          !klev

  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   DO k=1, klev-1
   DO i=1, klon
    sed_dust3D(i,k)=sed_dust3D(i,k) + &
          sed_flux(i,k+1)/zdz(i,k)
   ENDDO          !klon
   ENDDO          !klev
  !
  DO k = 1, klev
  DO i = 1, klon
     tr_seri(i,k,id_codu)=tr_seri(i,k,id_codu)+ &
           sed_dust3D(i,k)*time_step
  ENDDO
  ENDDO


   DO i=1, klon
     sed_dust(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
   ENDDO          !klon
   ELSE
    DO i=1, klon
      sed_dust(i)=0.
    ENDDO
   ENDIF
  !


  !--------- For scoarse  dust ------------------
  !
  !
   IF(id_scdu>0) THEN
   DO k=1, klev
   DO i=1,klon
  !
    v_stokes=RG*(rho_dust-zrho(i,k))* & !m/sec
          mmd_dustsco*mmd_dustsco* &
          1.e-12/(18.0*air_visco(i,k)/10.)
  !
   lambda=6.6*1.e-8*(103125/pplay(i,k))*(t_seri(i,k)/293.15)
   CC=1.0+1.257*lambda/(mmd_dustsco)/1.e6        !dimensionless
   v_sed=v_stokes*CC                       !m/sec
  !
  !---------check for v_sed*dt<zdz


   IF (v_sed*time_step.GT.zdz(i,k)) THEN
     v_sed=zdz(i,k)/time_step
   ENDIF

  !
   v_dep_dustsco(i,k)= v_sed
   sed_flux(i,k)     = tr_seri(i,k,id_scdu)*v_sed !g/cm3.m/sec
   ! !sed_dustsco3D(i,k)= -sed_flux(i,k)/zdz(i,k)      !g/cm3*sec !!!!!!!
  !
   ENDDO          !klon
   ENDDO          !klev

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   sed_dustsco3D(:,:)=0.0  ! initialisation

   DO k=1, klev
   DO i=1, klon
   sed_dustsco3D(i,k)=sed_dustsco3D(i,k)- &
         sed_flux(i,k)/zdz(i,k)
   ENDDO          !klon
   ENDDO          !klev
  !
   DO k=1, klev-1
   DO i=1, klon
    sed_dustsco3D(i,k)=sed_dustsco3D(i,k) + &
          sed_flux(i,k+1)/zdz(i,k)
   ENDDO          !klon
   ENDDO          !klev

  DO k = 1, klev
  DO i = 1, klon
   tr_seri(i,k,id_scdu)=tr_seri(i,k,id_scdu)+ &
         sed_dustsco3D(i,k)*time_step
  ENDDO
  ENDDO
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  !
   DO i=1, klon
     sed_dustsco(i)=sed_flux(i,1)*1.e6*1.e3    !--unit mg/m2/s
   ENDDO          !klon
   ELSE
    DO i=1, klon
      sed_dustsco(i)=0.
    ENDDO
   ENDIF
  !




  !
   RETURN
END SUBROUTINE sediment_mod
