subroutine deposition_source(ngrid, nlayer, nq,   &
                           ig, zzlay, zzlev,zdens,      &
                           zycol, ptimestep)
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!     dry deposition of chemical species
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
      use chimiedata_h
      use gases_h
      use tracer_h, only: noms, nesp, mmol
      use conc_mod, only: mmean ! mean molecular mass of the atmosphere

      implicit none
!
!
!     input
! 
      integer,intent(in) :: ngrid    ! number of atmospheric columns
      integer,intent(in) :: nlayer   ! number of atmospheric layers
      integer,intent(in) :: nq       ! number of tracers
      integer ig                     ! grid point index
      real    zzlay(ngrid,nlayer)    ! altitude at the middle of the layers (m)
      real    zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries (m)
      real    zdens(nlayer)          ! density (cm^-3)
      real    zycol(nlayer,nesp)     ! composition (volume mixing ratio)
      real    ptimestep              ! physical timestep (s)
!
!     local
!
      real    vd                         ! dry deposition velocity (cm.s-1)
      real    deltaz                     ! thickness of first layer (m)
      real    deltaznlay                 ! thickness of last layer (m)
      real    loss                       ! loss rate (s-1)
      real    prod                       ! production rate (s-1) 
      integer iq
      logical, save :: firstcall = .true.


      if (firstcall) then
         print*,'photochemistry: initialize deposition/source'

         ! You can set SF_mode/SF_value/prod_rate in traceur.def with #Moderntrac-v1
         !! Cases SF_mode=1 (fixed mixing ratio)
         !SF_mode(igcm_co2)=1
         !SF_value(igcm_co2)=gfrac(igas_CO2)

         !! Cases SF_mode=2 (production flux in molecules/m2/s)
         !prod_rate(indexchim('co'))=3.0e13

         firstcall=.false.
      endif !firstcall

      ! thickness of first layer (m)
      deltaz = zzlev(ig,2) - zzlev(ig,1)
      deltaznlay = zzlev(ig,nlayer) - zzlev(ig,nlayer-1)

      ! hydrogen escape
      !escape(ig,indexchim('h2'))=2.5e17*(zycol(nlayer,indexchim('h2'))+2*zycol(nlayer,indexchim('ch4'))+zycol(nlayer,indexchim('h2o_vap')))

      do iq=1,nesp
        if(SF_mode(iq).eq.1) then
          surface_flux(ig,iq)=(SF_value(iq)-zycol(1,iq))/ptimestep*(zdens(1)*1e6*deltaz)
          zycol(1,iq) = SF_value(iq)
          !! Set mmr instead of vmr
          !surface_flux(ig,iq)=(SF_value(iq)*mmean(ig,1)/mmol(iq+1)-zycol(1,iq))/ptimestep*(zdens(1)*1e6*deltaz)
          !zycol(1,iq) = SF_value(iq)*mmean(ig,1)/mmol(iq+1)
          zycol(nlayer,iq) = zycol(nlayer,iq)-ptimestep*escape(ig,iq)/(zdens(nlayer)*1e6*deltaznlay)
        elseif(SF_mode(iq).eq.2) then ! loss/prod rate (s-1)
          loss = 0.01*SF_value(iq)/deltaz
          zycol(1,iq) = zycol(1,iq)*exp(-loss*ptimestep)+(prod_rate(iq)/(zdens(1)*1e6*deltaz))*ptimestep
          zycol(nlayer,iq) = zycol(nlayer,iq)-ptimestep*escape(ig,iq)/(zdens(nlayer)*1e6*deltaznlay)
          surface_flux(ig,iq)=-0.01*SF_value(iq)*zycol(1,iq)*(zdens(1)*1e6)+prod_rate(iq)
        endif
      enddo ! end nesp

      return
      end