MODULE module_dry_dep


  USE chem1_list, ONLY:           &
       nspecies_chem =>nspecies,  & ! (IN)
       spc_alloc_chem=>spc_alloc, & ! (IN)
       spc_name_chem=>spc_name,   & ! (IN)
       ddp,                       & ! (IN)
       transport,                 & ! (IN)
       off,                       & ! (IN)
!DSM       O3,                        & ! (IN)
       dvj,                       & ! (IN)
       hstar,                     & ! (IN)
       ak0,                       & ! (IN)
       dak,                       & ! (IN)
       dhr,                       & ! (IN)
       f0                           ! (IN)

  USE aer1_list, ONLY:            &
       nspecies_aer=> nspecies,   & ! (IN)
       spc_alloc_aer=>spc_alloc,  & ! (IN)
       part_radius,               & ! (IN)
       part_dens,                 & ! (IN)
       mode_alloc,                & ! (IN)
       on,                        & ! (IN)
       nmodes                       ! (IN)

  USE mem_chem1, ONLY:            &
       chem1_vars,                & ! Type
       chem1_g                      ! (INOUT)

  USE mem_aer1, ONLY:             &
       aer1_vars,                 & ! Type
       aer1_g                       ! (INOUT)

!--(DMK-CCATT-BRAMS-5.0-INI)------------------------------------------------------------------
  USE ModDateUtils, ONLY:         &
       julday
!--(DMK-CCATT-BRAMS-5.0-FIM)------------------------------------------------------------------


  IMPLICIT NONE

  PRIVATE


  INTEGER, PARAMETER :: dep_seasons = 5 
  INTEGER, PARAMETER :: nlu = 25
  INTEGER, PARAMETER :: nseason = 1
  INTEGER, PARAMETER :: nseasons = 2

  ! following currently hardwired to leaf-3 classes
  INTEGER, PARAMETER :: isice=2
  INTEGER, PARAMETER :: iswater=1
  !- aerosol section
  INTEGER, PARAMETER :: ntotal=nspecies_aer*nmodes
  INTEGER, PUBLIC :: naer_a,naer_z
  INTEGER, PUBLIC :: ind_aer(ntotal)
  INTEGER, PUBLIC :: ind_mode(ntotal)
  INTEGER, PUBLIC :: NAER_TRANSPORTED
  INTEGER, PUBLIC :: NAER_TRANSPORTED_numb
  
  TYPE, PUBLIC :: sedim_type

     REAL,POINTER :: v_sed_part(:,:,:)
     REAL,POINTER :: r_lsl_part(:,:,:)
     REAL,POINTER :: v_dep_part(:,:)
       
  END TYPE sedim_type

  TYPE(sedim_type), PUBLIC, ALLOCATABLE :: dd_sedim(:,:)
  TYPE(sedim_type), PUBLIC, ALLOCATABLE :: dd_sedim_numb(:)
  LOGICAL :: aer_alloc = .FALSE., aer_alloc_numb = .FALSE.
  
  INTEGER :: ixxxlu(nlu)
  REAL    :: kpart(nlu)
  REAL    :: rac(nlu,dep_seasons)
  REAL    :: rclo(nlu,dep_seasons)
  REAL    :: rcls(nlu,dep_seasons)
  REAL    :: rgso(nlu,dep_seasons)
  REAL    :: rgss(nlu,dep_seasons)
  REAL    :: ri(nlu,dep_seasons)
  REAL    :: rlu(nlu,dep_seasons)
  REAL    :: dratio(nspecies_chem)

  REAL,    EXTERNAL :: rs        ! Function

  PUBLIC :: dep_init             ! Subroutine
  PUBLIC :: alloc_aer_sedim
  PUBLIC :: alloc_aer_sedim_numb
  PUBLIC :: fa_preptc_with_sedim ! Subroutine
  PUBLIC :: dry_dep              ! Subroutine


CONTAINS

  !=========================================================================
  SUBROUTINE alloc_aer_sedim(npatch,ngrids &
                            ,nmodes,nspecies_aer    &
                            ,mode_alloc,on          &
                            ,mmzp,mmxp,mmyp)

    INTEGER , INTENT(IN) :: npatch
    INTEGER , INTENT(IN) :: ngrids

    ! aer1_list
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: nspecies_aer
    INTEGER , INTENT(IN) :: mode_alloc(nmodes,nspecies_aer)
    INTEGER , INTENT(IN) :: on

    ! node_mod
    INTEGER , INTENT(IN) :: mmxp(ngrids)
    INTEGER , INTENT(IN) :: mmyp(ngrids)
    INTEGER , INTENT(IN) :: mmzp(ngrids)

    INTEGER :: ispc,imode,ng


    IF(aer_alloc) THEN
       PRINT *,'ERROR: aer_alloc already allocated'
       PRINT *,'Routine: aer_alloc File: chem_dry_dep.f90'
       STOP
    END IF

    NAER_TRANSPORTED = 0

    ! aerosol section : mapping  
    naer_a = 1
    naer_z = 0

    DO ispc = 1,nspecies_aer
       DO imode=1,nmodes

          IF(mode_alloc   (imode,ispc) == on ) THEN 

             naer_z = naer_z + 1
             ind_aer (naer_z) = ispc
             ind_mode(naer_z) = imode

          ENDIF
       ENDDO
    ENDDO

    ! total number of species (aer) to be transported 
    NAER_TRANSPORTED  =  naer_z
    IF(NAER_TRANSPORTED == 0) RETURN

    ALLOCATE (dd_sedim(NAER_TRANSPORTED,ngrids))
    DO ng=1,ngrids
       DO ispc=1,NAER_TRANSPORTED
          ALLOCATE(dd_sedim(ispc,ng)%v_sed_part (mmzp(ng),mmxp(ng),mmyp(ng)))
          dd_sedim(ispc,ng)%v_sed_part = 0.
          ALLOCATE(dd_sedim(ispc,ng)%r_lsl_part (mmxp(ng),mmyp(ng),npatch))
          dd_sedim(ispc,ng)%r_lsl_part = 0.
          ALLOCATE(dd_sedim(ispc,ng)%v_dep_part (mmxp(ng),mmyp(ng)))
          dd_sedim(ispc,ng)%v_dep_part = 0.
       ENDDO
    ENDDO


    aer_alloc=.TRUE.

  END SUBROUTINE alloc_aer_sedim

  !=========================================================================
  SUBROUTINE alloc_aer_sedim_numb(npatch,ngrids &
                            ,nmodes  &
                            ,mode_alloc,on          &
                            ,mmzp,mmxp,mmyp)

    INTEGER , INTENT(IN) :: npatch
    INTEGER , INTENT(IN) :: ngrids

    ! aer1_list
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: mode_alloc(nmodes)
    INTEGER , INTENT(IN) :: on

    ! node_mod
    INTEGER , INTENT(IN) :: mmxp(ngrids)
    INTEGER , INTENT(IN) :: mmyp(ngrids)
    INTEGER , INTENT(IN) :: mmzp(ngrids)

    INTEGER :: ispc,imode,ng


    IF(aer_alloc_numb) THEN
       PRINT *,'ERROR: aer_alloc for number already allocated'
       PRINT *,'Routine: aer_alloc File: chem_dry_dep.f90'
       STOP
    END IF

    NAER_TRANSPORTED = 0

    ! aerosol section : mapping
    naer_a = 1
    naer_z = 0

    !DO ispc = 1,nspecies_aer
       DO imode=1,nmodes

          IF(mode_alloc   (imode) == on ) THEN

             naer_z = naer_z + 1
             ind_aer (naer_z) = ispc
             ind_mode(naer_z) = imode

          ENDIF
       ENDDO
    !ENDDO

    ! total number of species (aer) to be transported
    NAER_TRANSPORTED_numb  =  naer_z
    IF(NAER_TRANSPORTED_numb == 0) RETURN

    ALLOCATE (dd_sedim_numb(ngrids))
    DO ng=1,ngrids
       DO ispc=1,NAER_TRANSPORTED
          ALLOCATE(dd_sedim_numb(ng)%v_sed_part (mmzp(ng),mmxp(ng),mmyp(ng)))
          dd_sedim_numb(ng)%v_sed_part = 0.
          ALLOCATE(dd_sedim_numb(ng)%r_lsl_part (mmxp(ng),mmyp(ng),npatch))
          dd_sedim_numb(ng)%r_lsl_part = 0.
          ALLOCATE(dd_sedim_numb(ng)%v_dep_part (mmxp(ng),mmyp(ng)))
          dd_sedim_numb(ng)%v_dep_part = 0.
       ENDDO
    ENDDO


    aer_alloc_numb=.TRUE.

  END SUBROUTINE alloc_aer_sedim_numb

  !========================================================================
  SUBROUTINE dry_dep(m1,m2,m3,ia,iz,ja,jz                 &
                    ,cpi,cpor,p00,g,vonk                  &
                    ,jdim,dzt,zt,nzpmax                   &
                    ,npatch,dt,level,nnqparm              &
                    ,imonth1,idate1,iyear1                &
                    ,chemistry,aerosol                    &
                    ,theta,rv,pp,dn0,pi0,up,vp,tke        &
                    ,sflux_t,sflux_r,sflux_u,sflux_v      &  
                    ,r_aer,ustar,tstar,patch_area,veg,Z0m &
                    ,rcp,pcpg,rtgt,rshort,conprr          &
                    ,aer1_g                               &
                    ,chem1_g                              &
                    ,dd_sedim)

    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3
    INTEGER , INTENT(IN)    :: ia
    INTEGER , INTENT(IN)    :: iz
    INTEGER , INTENT(IN)    :: ja
    INTEGER , INTENT(IN)    :: jz

    ! rconstants
    REAL    , INTENT(IN)    :: cpi
    REAL    , INTENT(IN)    :: cpor
    REAL    , INTENT(IN)    :: p00
    REAL    , INTENT(IN)    :: g
    REAL    , INTENT(IN)    :: vonk

    ! mem_grid
    INTEGER , INTENT(IN)    :: npatch
    INTEGER , INTENT(IN)    :: jdim
    REAL    , INTENT(IN)    :: dt
    REAL    , INTENT(IN)    :: dzt(nzpmax)
    REAL    , INTENT(IN)    :: zt(nzpmax)
    INTEGER , INTENT(IN)    :: nzpmax
    INTEGER , INTENT(IN)    :: imonth1
    INTEGER , INTENT(IN)    :: idate1
    INTEGER , INTENT(IN)    :: iyear1

    ! micphys
    INTEGER , INTENT(IN)    :: level

    ! mem_basic
    REAL    , INTENT(IN)    :: theta(m1,m2,m3)
    REAL    , INTENT(IN)    :: rv(m1,m2,m3)
    REAL    , INTENT(IN)    :: pp(m1,m2,m3)
    REAL    , INTENT(IN)    :: dn0(m1,m2,m3)
    REAL    , INTENT(IN)    :: pi0(m1,m2,m3)
    REAL    , INTENT(IN)    :: up(m1,m2,m3)
    REAL    , INTENT(IN)    :: vp(m1,m2,m3)

    ! mem_turb
    REAL    , INTENT(IN)    :: tke(m1,m2,m3)
    REAL    , INTENT(IN)    :: sflux_t(m2,m3)
    REAL    , INTENT(IN)    :: sflux_r(m2,m3)
    REAL    , INTENT(IN)    :: sflux_u(m2,m3)
    REAL    , INTENT(IN)    :: sflux_v(m2,m3)

    ! mem_radiate
    REAL    , INTENT(IN)    :: rshort(m2,m3)

    ! mem_micro
    REAL    , INTENT(IN)    :: rcp(m1,m2,m3)
    REAL    , INTENT(IN)    :: pcpg(m2,m3)


    ! mem_grid
    REAL    , INTENT(IN)    :: rtgt(m2,m3)

    ! mem_leaf
    REAL    , INTENT(IN)    :: r_aer(m2,m3,npatch)
    REAL    , INTENT(IN)    :: ustar(m2,m3,npatch)
    REAL    , INTENT(IN)    :: tstar(m2,m3,npatch)
    REAL    , INTENT(IN)    :: patch_area(m2,m3,npatch)
    REAL    , INTENT(IN)    :: veg(m2,m3,npatch)
    REAL    , INTENT(IN)    :: Z0m(m2,m3,npatch)

    ! mem_cuparm
    INTEGER , INTENT(IN)    :: nnqparm
    REAL    , INTENT(IN)    :: conprr(m2,m3)

    ! mem_aer1
    INTEGER         , INTENT(IN)    :: aerosol
    TYPE(aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,nspecies_aer)

    ! mem_chem1
    INTEGER         , INTENT(IN)    :: chemistry
    TYPE(chem1_vars), INTENT(INOUT) :: chem1_g(nspecies_chem)

    TYPE(sedim_type), INTENT(INOUT) :: dd_sedim(naer_transported)


    REAL, PARAMETER ::  ubmin = 0.25

    INTEGER :: i,j,ipatch,idry_part,ispc,k
    REAL    :: ups,vps,pis,sbf

    REAL    :: v_dep_gas(nspecies_chem,m2,m3)
    REAL    :: check_rain(m2,m3              )
    REAL    :: rmol(m2,m3              )!1./Monin-Obukhob
    REAL    :: rhchem(m2,m3)     

    ! (DMK) old scratches vars
    REAL    :: temps(m2,m3), prss(m2,m3), dens(m2,m3), vels(m2,m3), rvs(m2,m3), Zi(m2,m3)
    REAL    :: temp3d(m1,m2,m3), air_dens3d(m1,m2,m3)

    !- no tracers, aerosols or chemistry, go back
    IF(CHEMISTRY < 0) RETURN

    v_dep_gas = 0.0   
    check_rain = 0.
    rmol = 0.
    !-aux variables

    DO j = ja,jz
       DO i = ia,iz
          rvs  (i,j) = rv(2,i,j)
          pis        = ( pp(1,i,j) + pp(2,i,j) + pi0(1,i,j) + pi0(2,i,j) ) * .5 * cpi  
          prss (i,j) = pis ** cpor * p00                                               
          dens (i,j) = ( dn0(1,i,j) + dn0(2,i,j) ) * .5
          temps(i,j) = theta(2,i,j) * pis        ! temps=theta*Exner/CP
          rhchem(i,j)=100.*MIN(1.,MAX(0.05,rvs(i,j)/rs(prss(i,j),temps(i,j))))
          ups        = (up(2,i,j) + up(2,i-1,j   )) * .5
          vps        = (vp(2,i,j) + vp(2,i,j-jdim)) * .5
          vels (i,j) = MAX(ubmin,SQRT(ups** 2 + vps** 2))
          !- compute surface buoyancy flux [sbf] and 1 / Monin Obukhov height [1/zl]	    
          sbf = (g/theta(2,i,j) * (1. + .61 * rv(2,i,j))) &
               * ( sflux_t(i,j) * (1. + .61 * rv(2,i,j)) &
               +   sflux_r(i,j) * .61 * theta(2,i,j)   )
          !-  reciprocal of the Monin-Obukhov length (1/m)
          rmol (i,j) = - vonk*sbf / MAX(0.1,SQRT(SQRT(sflux_u(i,j)**2 + sflux_v(i,j)**2)))
          IF(ABS(rmol (i,j)) < 1.e-7) rmol (i,j) = 0. 
          !    print*,'rmol=',rmol(i,j)
       ENDDO
    ENDDO
    IF (nnqparm > 0) THEN
       check_rain(ia:iz,ja:jz)  = conprr(ia:iz,ja:jz)    
    ENDIF
    IF (level >= 3) THEN
       check_rain(ia:iz,ja:jz)  = check_rain(ia:iz,ja:jz)  +  pcpg(ia:iz,ja:jz)
    ENDIF

    CALL define_PBL_height(m1,m2,m3,npatch,ia,iz,ja,jz,zt     &
         ,tke,sflux_t,rcp,rtgt &
         ,zi)


    !---  dry deposition for gases
!DSM    CALL dry_dep_gases(m1,m2,m3,nspecies_chem,npatch,ia,iz,ja,jz  &
!DSM         ,v_dep_gas,r_aer  		           &
!DSM         ,prss,temps,dens,vels,rvs,rcp,Zi		   &
!DSM         ,ustar,tstar,patch_area,veg,Z0m		   &
!DSM         ,rshort,rtgt,dzt,check_rain,rmol,rhchem &
!DSM         ,O3,spc_alloc_chem,transport &
!DSM         ,on,dvj,hstar,ak0,dak,dhr,f0,imonth1,idate1,iyear1)
    CALL dry_dep_gases(m1,m2,m3,nspecies_chem,npatch,ia,iz,ja,jz  &
         ,v_dep_gas,r_aer  		           &
         ,prss,temps,dens,vels,rvs,rcp,Zi		   &
         ,ustar,tstar,patch_area,veg,Z0m		   &
         ,rshort,rtgt,dzt,check_rain,rmol,rhchem &
         ,spc_alloc_chem,transport &
         ,on,dvj,hstar,ak0,dak,dhr,f0,imonth1,idate1,iyear1)

    !-apply dry deposition on the tracers concentration and get the deposited mass on surface
    DO ispc=1,nspecies_chem
       IF (spc_alloc_chem(ddp,ispc) == off .OR. &
           spc_alloc_chem(transport,ispc) == off) CYCLE

       CALL apply_drydep(m1,m2,m3,ia,iz,ja,jz		&
            ,v_dep_gas(ispc,1:m2,1:m3)	& ! don't change this
            ,chem1_g  (ispc)%sc_t     &! tendency array
            ,chem1_g  (ispc)%sc_dd    &! deposited mass
            ,chem1_g  (ispc)%sc_p     &! mixing ratio
            ,dens,rtgt,dzt,dt)       
    ENDDO

    !--- aerosol dry deposition and sedimentation 

    IF(AEROSOL > 0) THEN

       DO j = ja,jz
          DO i = ia,iz
             DO k = 1, m1-1
                temp3d(k,i,j) = 0.5*( theta(k+1,i,j) * ( pp(k+1,i,j)+pi0(k+1,i,j) ) * cpi + &
                     theta(k  ,i,j) * ( pp(k  ,i,j)+pi0(k  ,i,j) ) * cpi   )
                air_dens3d(k,i,j) = 0.5 *(  dn0(k+1,i,j) +  dn0(k  ,i,j) )
             ENDDO
             temp3d(m1,i,j) =	 temp3d(m1-1,i,j) 

             air_dens3d(m1,i,j) =    air_dens3d(m1-1,i,j)
          ENDDO
       ENDDO
       !print*,'temp=',temp3d(1:2,int(iz/2),int(jz/2))

       !--  dry deposition/sedimentation for particles
       CALL dry_dep_sedim_particles(m1,m2,m3,nspecies_aer              &
                                   ,npatch,ia,iz,ja,jz                 &
                                   ,r_aer,temp3d,air_dens3d,temps,dens &
                                   ,vels,rvs,Zi,ustar,tstar,patch_area &
                                   ,veg,Z0m,nmodes,nspecies_aer        &
                                   ,part_radius,part_dens,mode_alloc   &
                                   ,on,g,dd_sedim)

       !- loop over the transported aerosols particles
       DO ispc=naer_a,naer_z


          IF(spc_alloc_aer(ddp      ,ind_mode(ispc),ind_aer(ispc)) == off .OR. &
               spc_alloc_aer(transport,ind_mode(ispc),ind_aer(ispc)) == off) CYCLE


          !print*,'dd=',dd_sedim(ispc,ngrid)%v_dep_part(int(iz/2),int(jz/2))*100.

          CALL apply_drydep(m1,m2,m3,ia,iz,ja,jz	       &
               ,dd_sedim(ispc)%v_dep_part   	               &
               ,aer1_g  (ind_mode(ispc),ind_aer(ispc))%sc_t    &! tendency array
               ,aer1_g  (ind_mode(ispc),ind_aer(ispc))%sc_dd   &! deposited mass
               ,aer1_g  (ind_mode(ispc),ind_aer(ispc))%sc_p    &! mixing ratio
               ,dens,rtgt,dzt,dt)

       ENDDO
    ENDIF  ! endif of aerosol

  END SUBROUTINE dry_dep
  !========================================================================

!DSM  SUBROUTINE dry_dep_gases(m1,m2,m3,nspecies_chem,npatch,ia,iz,ja,jz         &
!DSM                          ,V_dep,r_aer,prss,temps,dens,vels,rvs,rcp,Zi       &
!DSM                          ,ustar,tstar,patch_area,veg,Z0m,rshort,rtgt,dzt    &
!DSM                          ,check_rain,rmol,rhchem,O3,spc_alloc_chem,transport&
!DSM                          ,on,dvj,hstar,ak0,dak,dhr,f0,imonth1,idate1,iyear1)
  SUBROUTINE dry_dep_gases(m1,m2,m3,nspecies_chem,npatch,ia,iz,ja,jz         &
                          ,V_dep,r_aer,prss,temps,dens,vels,rvs,rcp,Zi       &
                          ,ustar,tstar,patch_area,veg,Z0m,rshort,rtgt,dzt    &
                          ,check_rain,rmol,rhchem,spc_alloc_chem,transport&
                          ,on,dvj,hstar,ak0,dak,dhr,f0,imonth1,idate1,iyear1)

    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3
    INTEGER , INTENT(IN)    :: nspecies_chem
    INTEGER , INTENT(IN)    :: npatch
    INTEGER , INTENT(IN)    :: ia
    INTEGER , INTENT(IN)    :: iz
    INTEGER , INTENT(IN)    :: ja
    INTEGER , INTENT(IN)    :: jz
    REAL    , INTENT(INOUT) :: V_dep(nspecies_chem,m2,m3)
    REAL    , INTENT(IN)    :: r_aer(m2,m3,npatch)
    REAL    , INTENT(IN)    :: prss(m2,m3)
    REAL    , INTENT(IN)    :: temps(m2,m3)

    REAL    , INTENT(IN)    :: dens(m2,m3)         ! (DMK) not used
    REAL    , INTENT(IN)    :: vels(m2,m3)         ! (DMK) not used
    REAL    , INTENT(IN)    :: rvs(m2,m3)          ! (DMK) not used
    REAL    , INTENT(IN)    :: rcp(m1,m2,m3)
    REAL    , INTENT(IN)    :: Zi(m2,m3)           ! (DMK) not used
    REAL    , INTENT(IN)    :: ustar(m2,m3,npatch)
    REAL    , INTENT(IN)    :: tstar(m2,m3,npatch) ! (DMK) not used
    REAL    , INTENT(IN)    :: patch_area(m2,m3,npatch)
    REAL    , INTENT(IN)    :: veg(m2,m3,npatch)
    REAL    , INTENT(IN)    :: Z0m(m2,m3,npatch)   ! (DMK) not used
    REAL    , INTENT(IN)    :: rshort(m2,m3)
    REAL    , INTENT(IN)    :: rtgt(m2,m3)
    REAL    , INTENT(IN)    :: dzt(m1)
    REAL    , INTENT(IN)    :: check_rain(m2,m3)
    REAL    , INTENT(IN)    :: rmol(m2,m3)
    REAL    , INTENT(IN)    :: rhchem(m2,m3)

    ! chem1_list
!DSM    INTEGER , INTENT(IN)    :: O3
    INTEGER , INTENT(IN)    :: spc_alloc_chem(6,nspecies_chem)
    INTEGER , INTENT(IN)    :: transport
    INTEGER , INTENT(IN)    :: on
    REAL    , INTENT(IN)    :: dvj(nspecies_chem)
    REAL    , INTENT(IN)    :: hstar(nspecies_chem)
    REAL    , INTENT(IN)    :: ak0(nspecies_chem)
    REAL    , INTENT(IN)    :: dak(nspecies_chem)
    REAL    , INTENT(IN)    :: dhr(nspecies_chem)
    REAL    , INTENT(IN)    :: f0(nspecies_chem)
    
    ! mem_grid
    INTEGER , INTENT(IN)    :: imonth1
    INTEGER , INTENT(IN)    :: idate1
    INTEGER , INTENT(IN)    :: iyear1

    INTEGER :: i,j,ipatch

    !- local var
    REAL :: vgs0d(nspecies_chem)  !,srfres
    INTEGER iveg,iseason,jday
    LOGICAL :: highnh3, rainflag, vegflag, wetflag
    REAL dvpart,dvfog,clwchem,z1,ustarw,vgpart
    !real , dimension(m2,m3) :: aer_res
    !real   ::  r_lsl(m2,m3,npatch,nspecies_chem),rcx(npatch,nspecies_chem)
    REAL   ::  r_lsl(nspecies_chem),rcx(nspecies_chem)

    !   Set the reference height (10.0 m)
    REAL, PARAMETER ::  zr = 10.0


    jday=julday(imonth1,idate1,iyear1)
    iseason = 1
    IF( jday.LT.90 .OR. jday.GT.270)  iseason=2

    DO j = ja,jz
       DO i = ia,iz
          ustarw = 0.
          !
          !     Set logical default values
          rainflag = .FALSE.
          IF (check_rain(i,j) > 0.  ) rainflag = .TRUE.

          wetflag  = .FALSE.
          IF (rhchem(i,j) >= 95.) wetflag  = .TRUE.

          clwchem = rcp(2,i,j)

          !-not using
          highnh3 = .FALSE.
          !      if(chem(i,kts,j,p_nh3).gt.2.*chem(i,kts,j,p_so2))highnh3 = .true.
          !----


          !zntt = znt(i,j)  ! Surface roughness height (m)

          z1 =   rtgt(i,j)/dzt(2)! height of the fisrt cell

          DO ipatch = 1,npatch
             IF (patch_area(i,j,ipatch) .GE. .009) THEN 

                IF(ipatch == 1) THEN
                   iveg = 1           ! rc routine needs water as 1
                ELSE
                   iveg = NINT(veg(i,j,ipatch))     
                ENDIF

                ustarw=ustarw+ustar(i,j,ipatch)*patch_area(i,j,ipatch)	
                !
                !- calculates surface resistance (rcx) 
                !

                !DSM CALL rc(rcx,temps(i,j),rshort(i,j),rhchem(i,j),iveg,iseason   &
                !DSM      ,nspecies_chem,wetflag,rainflag,highnh3,spc_alloc_chem   &
                !DSM      ,hstar,ak0,dak,dhr,f0,transport,O3,on)

                CALL rc(rcx,temps(i,j),rshort(i,j),rhchem(i,j),iveg,iseason   &
                     ,nspecies_chem,wetflag,rainflag,highnh3,spc_alloc_chem   &
                     ,hstar,ak0,dak,dhr,f0,transport,on)

                !- get R_b 
                CALL get_rb(nspecies_chem,r_lsl,temps(i,j) &
                     ,prss(i,j),ustar(i,j,ipatch),dvj)
                !
                !- deposition velocity for each patch 
                vgs0d(1:nspecies_chem) = 1./( r_aer(i,j,ipatch     )  + &
                     r_lsl(1:nspecies_chem)  + &
                     rcx  (1:nspecies_chem)    )
                !			  
                !- special treatment for SULF
!!$	    if( spc_alloc_chem(transport,SULF) == ON) then
!!$             CALL deppart(rmol(i,j),ustar(i,j,ipatch),rhchem(i,j),clwchem,iveg,dvpart,dvfog)
!!$             CALL depvel(nspecies_chem,rmol(i,j),zr,Z0m(i,j,ipatch)   &
!!$	               ,ustar(i,j,ipatch),vgpart)
!!$             vgs0d(sulf)=1.0/((1.0/vgpart)+(1.0/dvpart))
!!$	    endif
                !- end SULF treatment
                !
                !
                !- effective deposition velocity for all patches
                V_dep(1:nspecies_chem,i,j) = V_dep(1:nspecies_chem,i,j) + &
                     patch_area(i,j,ipatch)*vgs0d(1:nspecies_chem)

             ENDIF
          ENDDO
          !	print*,'O3-1 cm/s=',V_dep(O3,i,j)*100; call flush(6)
          !- get the effective deposition velocity
          CALL cellvg(V_dep(1:nspecies_chem,i,j),ustarw,z1,zr,rmol(i,j),nspecies_chem)
          !	print*,'O3-2 cm/s=',V_dep(O3,i,j)*100; call flush(6)

       ENDDO
    ENDDO
    RETURN
  END SUBROUTINE dry_dep_gases
  !========================================================================

  SUBROUTINE get_rb(nspecies_chem,r_lsl,temps,prss,ustar,dvj)

    INTEGER , INTENT(IN)  :: nspecies_chem
    REAL    , INTENT(OUT) :: r_lsl(nspecies_chem)
    REAL    , INTENT(IN)  :: temps
    REAL    , INTENT(IN)  :: prss
    REAL    , INTENT(IN)  :: ustar

    ! chem1_list
    REAL    , INTENT(IN)  :: dvj(nspecies_chem)

    !local var
    REAL, PARAMETER :: i23=2./3., STDTEMP = 273.15, STDATMPA = 101325.0 
    INTEGER ispc
    REAL dvj_tmp,sc,scpr23,eta


    !-  kinematic viscosity of air at surface (cm^2/s)
    eta = -1.1555E-10*temps**3 + 9.5728E-07*temps**2 + 3.7604E-04*temps - 3.4484E-02  

    DO ispc=1,nspecies_chem

       !         dvj(ispc) = dvj(ispc)*(293.15/298.15)**1.75
       !         dvj_tmp = dvj(ispc)*(temps/(temps+5.))**1.75
       dvj_tmp = dvj(ispc)* (( temps / STDTEMP ) ** 1.75)&
            * (STDATMPA / prss ) ! cm^2/s
       !          dratio(ispc) = 0.242/dvj(ispc) 

       !         sc = 0.15/dvj_tmp ! Schmidt Number at 20�C   
       sc =  eta/dvj_tmp ! Schmidt Number at any temps 

       scpr23 = (sc/0.72)**(i23) ! (Schmidt # / Prandtl #)**
       r_lsl(ispc) = 5.*scpr23/ustar
    ENDDO
  END SUBROUTINE get_rb

  ! **********************************************************************
  ! **************************  SUBROUTINE RC  ***************************
  ! **********************************************************************
!DSM   SUBROUTINE rc(rcx,t,rad,rh,iland,iseason,nspecies_chem,        &
!DSM                 wetflag,rainflag,highnh3,spc_alloc_chem,hstar,   &
!DSM                 ak0,dak,dhr,f0,transport,O3,on )
  SUBROUTINE rc(rcx,t,rad,rh,iland,iseason,nspecies_chem,        &
                wetflag,rainflag,highnh3,spc_alloc_chem,hstar,   &
                ak0,dak,dhr,f0,transport,on )
    !     THIS SUBROUTINE CALCULATES SURFACE RESISTENCES ACCORDING
    !     TO THE MODEL OF
    !     M. L. WESELY,
    !     ATMOSPHERIC ENVIRONMENT 23 (1989), 1293-1304
    !     WITH SOME ADDITIONS ACCORDING TO
    !     J. W. ERISMAN, A. VAN PUL, AND P. WYERS,
    !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607
    !     WRITTEN BY  WINFRIED SEIDL, APRIL 1997
    !     MODYFIED BY WINFRIED SEIDL, MARCH 2000
    !                    FOR MM5 VERSION 3
    !----------------------------------------------------------------------

    REAL    , INTENT(IN)    :: t
    REAL    , INTENT(IN)    :: rad
    REAL    , INTENT(IN)    :: rh ! (DMK) not used
    INTEGER , INTENT(IN)    :: iland   ! (DMK) IN
    INTEGER , INTENT(IN)    :: iseason ! (DMK) IN
    INTEGER , INTENT(IN)    :: nspecies_chem
    LOGICAL , INTENT(IN)    :: wetflag
    LOGICAL , INTENT(IN)    :: rainflag
    LOGICAL , INTENT(IN)    :: highnh3 ! (DMK) not used
    REAL    , INTENT(INOUT) :: rcx(nspecies_chem)
    
    ! chem1_list
    INTEGER , INTENT(IN)    :: spc_alloc_chem(6,nspecies_chem)
    REAL    , INTENT(IN)    :: hstar(nspecies_chem)
    REAL    , INTENT(IN)    :: ak0(nspecies_chem)
    REAL    , INTENT(IN)    :: dak(nspecies_chem)
    REAL    , INTENT(IN)    :: dhr(nspecies_chem)
    REAL    , INTENT(IN)    :: f0(nspecies_chem)
    INTEGER , INTENT(IN)    :: transport
!DSM    INTEGER , INTENT(IN)    :: O3
    INTEGER , INTENT(IN)    :: on

!------------old--------------------------------------------
! .. Scalar Arguments ..
!        REAL :: rad, rh, t, temp_i
!        INTEGER :: iland, iseason, nspecies_chem
!        LOGICAL :: highnh3, rainflag, wetflag
!        real :: rcx(nspecies_chem)
!------------old--------------------------------------------

    REAL :: temp_i
    ! ..
    ! .. Array Arguments ..
    INTEGER :: iprt
    ! ..
    ! .. Local Scalars ..
    REAL :: rclx, rdc, resice, rgsx, rluo1, rluo2, rlux, rmx, rs, rsmx, &
         tc, rdtheta, z, hplus, corrh
    INTEGER :: n
    ! ..
    ! .. Local Arrays ..
    REAL :: hstary(nspecies_chem)
    ! ..
    ! .. Intrinsic Functions ..
    INTRINSIC exp
    !srf- declaring NH3, because it is not in the chem1-list for RACM
    INTEGER :: NH3 ! see return comand before the calculation of RC for NH3
    !srf---
    ! ..
    INTEGER, PARAMETER :: SO2=10000

    DO n = 1, nspecies_chem
       rcx(n) = 1.
    END DO

    tc = t - 273.15
    rdtheta = 0.
    temp_i = 1./t-1./298.

    z = 200./(rad+0.1)

!!!  HARDWIRE VALUES FOR TESTING
    !       z=0.4727409
    !       tc=22.76083
    !       t=tc+273.15
    !       rad = 412.8426
    !       rainflag=.false.
    !       wetflag=.false.

    IF ((tc<=0.) .OR. (tc>=40.)) THEN
       rs = 9999.
    ELSE

!--(DMK-CCATT-INI)----------------------------------------------------------
       rs = ri(ixxxlu(iland),iseason)*(1+z*z)*(400./(tc*(40.-tc)))
!--(DMK-CCATT-OLD)----------------------------------------------------------
!       rs = ri(iland,iseason)*(1+z*z)*(400./(tc*(40.-tc)))
!--(DMK-CCATT-FIM)----------------------------------------------------------
    END IF
    rdc = 100*(1.+1000./(rad+10))/(1+1000.*rdtheta)

!--(DMK-CCATT-INI)----------------------------------------------------------
    rluo1 = 1./(1./3000.+1./3./rlu(ixxxlu(iland),iseason))
    rluo2 = 1./(1./1000.+1./3./rlu(ixxxlu(iland),iseason))
!--(DMK-CCATT-OLD)----------------------------------------------------------
!    rluo1 = 1./(1./3000.+1./3./rlu(iland,iseason))
!    rluo2 = 1./(1./1000.+1./3./rlu(iland,iseason))
!--(DMK-CCATT-FIM)----------------------------------------------------------

    resice = 1000.*EXP(-tc-4.)

    ! change MP 11/12/07 taking into account the acid dissociation constant
    ! and assuming pH=7 in the vegetation
    !srf-opt  hstary(n) = hstar(n)*exp(dhr(n)*(1./t-1./298.))
    !                    *(1+ak0(n)*exp(dak(n)*(1./t-1/298.))/hplus)
    !          pH=7  hplus=10**(-pH) 
    hplus=1.E-7
    DO n = 1, nspecies_chem
       !srf      IF (hstar(n)==0.) GO TO 10
       IF (hstar(n)==0.) CYCLE
       !srf
       corrh=1+ak0(n)*EXP(dak(n)*(temp_i))/hplus
       hstary(n) = hstar(n)*EXP(dhr(n)*(temp_i))*corrh
       !
       !         if ((n == 8).or.(n == 52))then
       !           print*, n, hstary(n),t,temp_i,ak0(n), dak(n),hplus,corrh &
       !                , hstar(n), dhr(n)
       !         endif
       ! end change MP
       rmx = 1./(hstary(n)/3000.+100.*f0(n))
       rsmx = rs*dratio(n) + rmx

!--(DMK-CCATT-INI)----------------------------------------------------------
       rclx = 1./(hstary(n)/1.E+5/rcls(ixxxlu(iland),iseason)+f0(n)/rclo(ixxxlu(iland), &
            iseason)) + resice
!--(DMK-CCATT-OLD)----------------------------------------------------------
!       rclx = 1./(hstary(n)/1.E+5/rcls(iland,iseason)+f0(n)/rclo(iland, &
!            iseason)) + resice
!--(DMK-CCATT-FIM)----------------------------------------------------------

!--(DMK-CCATT-INI)----------------------------------------------------------
       rgsx = 1./(hstary(n)/1.E+5/rgss(ixxxlu(iland),iseason)+f0(n)/rgso(ixxxlu(iland), &
            iseason)) + resice
!--(DMK-CCATT-OLD)----------------------------------------------------------
!       rgsx = 1./(hstary(n)/1.E+5/rgss(iland,iseason)+f0(n)/rgso(iland, &
!            iseason)) + resice
!--(DMK-CCATT-FIM)----------------------------------------------------------

!--(DMK-CCATT-INI)----------------------------------------------------------
       rlux = rlu(ixxxlu(iland),iseason)/(1.E-5*hstary(n)+f0(n)) + resice
!--(DMK-CCATT-OLD)----------------------------------------------------------
!       rlux = rlu(iland,iseason)/(1.E-5*hstary(n)+f0(n)) + resice
!--(DMK-CCATT-FIM)----------------------------------------------------------

       IF (wetflag) THEN

!--(DMK-CCATT-INI)----------------------------------------------------------
          rlux = 1./(1./3./rlu(ixxxlu(iland),iseason)+1.E-7*hstary(n)+f0(n)/rluo1)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rlux = 1./(1./3./rlu(iland,iseason)+1.E-7*hstary(n)+f0(n)/rluo1)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END IF
       IF (rainflag) THEN

!--(DMK-CCATT-INI)----------------------------------------------------------
          rlux = 1./(1./3./rlu(ixxxlu(iland),iseason)+1.E-7*hstary(n)+f0(n)/rluo2)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rlux = 1./(1./3./rlu(iland,iseason)+1.E-7*hstary(n)+f0(n)/rluo2)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END IF

!--(DMK-CCATT-INI)----------------------------------------------------------
       rcx(n) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(ixxxlu(iland), &
            iseason)+rgsx))
!--(DMK-CCATT-OLD)----------------------------------------------------------
!       rcx(n) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(iland, &
!            iseason)+rgsx))
!--(DMK-CCATT-FIM)----------------------------------------------------------

       IF (rcx(n)<1.) rcx(n) = 1.
10  END DO
    !

!DSM    IF( spc_alloc_chem(transport,O3) == ON) THEN

       !--------------------------------- O3 
       !        SPECIAL TREATMENT FOR OZONE
       !srf       hstary(O3) = hstar(O3)*exp(dhr(O3)*(1./t-1./298.))
!DSM       hstary(O3) = hstar(O3)*EXP(dhr(O3)*(temp_i))
!DSM       rmx = 1./(hstary(O3)/3000.+100.*f0(O3))
!DSM       rsmx = rs*dratio(O3) + rmx

!DSM !--(DMK-CCATT-INI)----------------------------------------------------------
!DSM        rlux = rlu(ixxxlu(iland),iseason)/(1.E-5*hstary(O3)+f0(O3)) + resice
!DSM        rclx = rclo(ixxxlu(iland),iseason) + resice
!DSM        rgsx = rgso(ixxxlu(iland),iseason) + resice
!DSM !--(DMK-CCATT-OLD)----------------------------------------------------------
!DSM !       rlux = rlu(iland,iseason)/(1.E-5*hstary(O3)+f0(O3)) + resice
!DSM !       rclx = rclo(iland,iseason) + resice
!DSM !       rgsx = rgso(iland,iseason) + resice
!DSM !--(DMK-CCATT-FIM)----------------------------------------------------------

!DSM        IF (wetflag) rlux = rluo1
!DSM        IF (rainflag) rlux = rluo2

!DSM !--(DMK-CCATT-INI)----------------------------------------------------------
!DSM        rcx(O3) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(ixxxlu(iland), &
!DSM             iseason)+rgsx))
!DSM !--(DMK-CCATT-OLD)----------------------------------------------------------
!DSM !       rcx(O3) = 1./(1./rsmx+1./rlux+1./(rdc+rclx)+1./(rac(iland, &
!DSM !            iseason)+rgsx))
!DSM !--(DMK-CCATT-FIM)----------------------------------------------------------

!DSM       rcx(O3)=MAX(rcx(O3),1.) !- old: IF (rcx(O3)<1.) rcx(O3) = 1.
!DSM    ENDIF


!!$      if( spc_alloc_chem(transport,SO2) == ON) then
!!$
!!$!    	      SPECIAL TREATMENT FOR SO2 (Wesely)
!!$!    		HSTARY(SO2)=HSTAR(SO2)*EXP(DHR(SO2)*(1./T-1./298.))
!!$!    		RMX=1./(HSTARY(SO2)/3000.+100.*F0(SO2))
!!$!    		RSMX=RS*DRATIO(SO2)+RMX
!!$!    		RLUX=RLU(ILAND,ISEASON)/(1.E-5*HSTARY(SO2)+F0(SO2))
!!$!    	     &       +RESICE
!!$!    		RCLX=RCLS(ILAND,ISEASON)+RESICE
!!$!    		RGSX=RGSS(ILAND,ISEASON)+RESICE
!!$!    		IF ((wetflag).OR.(RAINFLAG)) THEN
!!$!    		  IF (ILAND.EQ.1) THEN
!!$!    		    RLUX=50.
!!$!    		  ELSE
!!$!    		    RLUX=100.
!!$!    		  END IF
!!$!    		END IF
!!$!    		RCX(SO2)=1./(1./RSMX+1./RLUX+1./(RDC+RCLX)
!!$!    	     &  	      +1./(RAC(ILAND,ISEASON)+RGSX))
!!$!    		IF (RCX(SO2).LT.1.) RCX(SO2)=1.
!!$!
!!$!
!!$!----------------------------------------- SO2 
!!$!    	      SO2 according to Erisman et al. 1994
!!$!    		R_STOM
!!$     		rsmx = rs*dratio(SO2)
!!$!    		R_EXT
!!$     		IF (tc>(-1.)) THEN
!!$     		  IF (rh<81.3) THEN
!!$     		    rlux = 25000.*exp(-0.0693*rh)
!!$     		  ELSE
!!$     		    rlux = 0.58E12*exp(-0.278*rh)
!!$     		  END IF
!!$     		END IF
!!$     		IF (((wetflag) .OR. (rainflag)) .AND. (tc>(-1.))) THEN
!!$     		  rlux = 1.
!!$     		END IF
!!$     		IF ((tc>=(-5.)) .AND. (tc<=(-1.))) THEN
!!$     		  rlux = 200.
!!$     		END IF
!!$     		IF (tc<(-5.)) THEN
!!$     		  rlux = 500.
!!$     		END IF
!!$!    		INSTEAD OF R_INC R_CL and R_DC of Wesely are used
!!$     		rclx = rcls(iland,iseason)
!!$!    		DRY SURFACE
!!$     		rgsx = 1000.
!!$!    		WET SURFACE
!!$     		IF ((wetflag) .OR. (rainflag)) THEN
!!$     		  IF (highnh3) THEN
!!$     		    rgsx = 0.
!!$     		  ELSE
!!$     		    rgsx = 500.
!!$     		  END IF
!!$     		END IF
!!$!    		WATER
!!$     		IF (iland==iswater) THEN
!!$     		  rgsx = 0.
!!$     		END IF
!!$!    		SNOW
!!$     		IF ((iseason==4) .OR. (iland==isice)) THEN
!!$     		  IF (tc>2.) THEN
!!$     		    rgsx = 0.
!!$     		  END IF
!!$     		  IF ((tc>=(-1.)) .AND. (tc<=2.)) THEN
!!$     		    rgsx = 70.*(2.-tc)
!!$     		  END IF
!!$     		  IF (tc<(-1.)) THEN
!!$     		    rgsx = 500.
!!$     		  END IF
!!$     		END IF
!!$!    		TOTAL SURFACE RESISTENCE
!!$     		IF ((iseason/=4) .AND. (ixxxlu(iland)/=1) .AND. (iland/=iswater) .AND. &
!!$     		    (iland/=isice)) THEN
!!$     		  rcx(SO2) = 1./(1./rsmx+1./rlux+1./(rclx+rdc+rgsx))
!!$     		ELSE
!!$     		  rcx(SO2) = rgsx
!!$     		END IF
!!$     		IF (rcx(SO2)<1.) rcx(SO2) = 1.
!!$!
!!$      endif


    !------------ NOT IN USE BELOW -----------------------------------------------
    !srf- dry dep for NH3 NOT in use :
    RETURN 
    !srf ---

    !     NH3 according to Erisman et al. 1994
    !       R_STOM
    rsmx = rs*dratio(NH3)
    !       GRASSLAND (PASTURE DURING GRAZING)
!kml    IF (ixxxlu(iland)==3) THEN
    IF (ixxxlu(iland)==2 .OR. ixxxlu(iland)==3 .OR. ixxxlu(iland)==4) THEN
       IF (iseason==1) THEN
          !           SUMMER
          rcx(NH3) = 1000.
       END IF
       IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
          !           WINTER, NO SNOW
          IF (tc>-1.) THEN
             IF (rad/=0.) THEN
                rcx(NH3) = 50.
             ELSE
                rcx(NH3) = 100.
             END IF
             IF ((wetflag) .OR. (rainflag)) THEN
                rcx(NH3) = 20.
             END IF
          END IF
          IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
             rcx(NH3) = 200.
          END IF
          IF (tc<(-5.)) THEN
             rcx(NH3) = 500.
          END IF
       END IF
    END IF
    !       AGRICULTURAL LAND (CROPS AND UNGRAZED PASTURE)
!kml    IF (ixxxlu(iland)==2) THEN
    IF (ixxxlu(iland)==5 .OR. ixxxlu(iland)==6) THEN
       IF (iseason==1) THEN
          !           SUMMER
          IF (rad/=0.) THEN
             rcx(NH3) = rsmx
          ELSE
             rcx(NH3) = 200.
          END IF
          IF ((wetflag) .OR. (rainflag)) THEN
             rcx(NH3) = 50.
          END IF
       END IF
       IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
          !           WINTER, NO SNOW
          IF (tc>-1.) THEN
             IF (rad/=0.) THEN
                rcx(NH3) = rsmx
             ELSE
                rcx(NH3) = 300.
             END IF
             IF ((wetflag) .OR. (rainflag)) THEN
                rcx(NH3) = 100.
             END IF
          END IF
          IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
             rcx(NH3) = 200.
          END IF
          IF (tc<(-5.)) THEN
             rcx(NH3) = 500.
          END IF
       END IF
    END IF
    !       SEMI-NATURAL ECOSYSTEMS AND FORESTS
!kml    IF ((ixxxlu(iland)==4) .OR. (ixxxlu(iland)==5) .OR. (ixxxlu( &
!kml         iland)==6)) THEN

    IF ((ixxxlu(iland)==10) .OR. (ixxxlu(iland)==11) .OR.  &
        (ixxxlu(iland)==12) .OR. (ixxxlu(iland)==13) .OR.  &
        (ixxxlu(iland)==14) .OR. (ixxxlu(iland)==15) .OR.  &
        (ixxxlu(iland)==18)) THEN
       IF (rad/=0.) THEN
          rcx(NH3) = 500.
       ELSE
          rcx(NH3) = 1000.
       END IF
       IF ((wetflag) .OR. (rainflag)) THEN
          IF (highnh3) THEN
             rcx(NH3) = 100.
          ELSE
             rcx(NH3) = 0.
          END IF
       END IF
       IF ((iseason==2) .OR. (iseason==3) .OR. (iseason==5)) THEN
          !           WINTER, NO SNOW
          IF ((tc>=(-5.)) .AND. (tc<=-1.)) THEN
             rcx(NH3) = 200.
          END IF
          IF (tc<(-5.)) THEN
             rcx(NH3) = 500.
          END IF
       END IF
    END IF
    !       WATER
    IF (iland==iswater) THEN
       rcx(NH3) = 0.
    END IF
    !       URBAN AND DESERT (SOIL SURFACES)
!kml    IF (ixxxlu(iland)==1) THEN
    IF (ixxxlu(iland)==1 .OR. ixxxlu(iland)==19 ) THEN
       IF ( .NOT. wetflag) THEN
          rcx(NH3) = 50.
       ELSE
          rcx(NH3) = 0.
       END IF
    END IF
    !       SNOW COVERED SURFACES OR PERMANENT ICE
    IF ((iseason==4) .OR. (iland==isice)) THEN
       IF (tc>2.) THEN
          rcx(NH3) = 0.
       END IF
       IF ((tc>=(-1.)) .AND. (tc<=2.)) THEN
          rcx(NH3) = 70.*(2.-tc)
       END IF
       IF (tc<(-1.)) THEN
          rcx(NH3) = 500.
       END IF
    END IF
    IF (rcx(NH3)<1.) rcx(NH3) = 1.

  END SUBROUTINE rc
  ! **********************************************************************
  SUBROUTINE cellvg(vgtemp,ustar,dz,zr,rmol,nspec)
    !     THIS PROGRAM HAS BEEN DESIGNED TO CALCULATE THE CELL AVERAGE
    !     DEPOSITION VELOCITY GIVEN THE VALUE OF VG AT SOME REFERENCE
    !     HEIGHT ZR WHICH IS MUCH SMALLER THAN THE CELL HEIGHT DZ.
    !       PROGRAM WRITTEN BY GREGORY J.MCRAE (NOVEMBER 1977)
    !         Modified by Darrell A. Winner    (February 1991)
    !.....PROGRAM VARIABLES...
    !     VgTemp   - DEPOSITION VELOCITY AT THE REFERENCE HEIGHT
    !     USTAR    - FRICTION VELOCITY
    !     RMOL     - RECIPROCAL OF THE MONIN-OBUKHOV LENGTH
    !     ZR       - REFERENCE HEIGHT
    !     DZ       - CELL HEIGHT
    !     CELLVG   - CELL AVERAGE DEPOSITION VELOCITY
    !     VK       - VON KARMAN CONSTANT

    REAL    , INTENT(IN)  :: ustar
    REAL    , INTENT(IN)  :: dz
    REAL    , INTENT(IN)  :: zr
    REAL    , INTENT(IN)  :: rmol
    INTEGER , INTENT(IN)  :: nspec
    REAL    , INTENT(OUT) :: vgtemp(nspec)

    ! ..
    ! .. Local Scalars ..
    REAL :: a, fac, pdz, pzr,fx
    INTEGER :: nss
    ! ..
    ! .. Intrinsic Functions ..
    INTRINSIC alog, sqrt
    ! ..
    !     Set the von Karman constant
    REAL, PARAMETER :: vk = 0.4

    !     DETERMINE THE STABILITY BASED ON THE CONDITIONS
    !             1/L < 0 UNSTABLE
    !             1/L = 0 NEUTRAL
    !             1/L > 0 STABLE


    IF (rmol<0) THEN
       pdz = SQRT(1.0-9.0*dz*rmol)
       pzr = SQRT(1.0-9.0*zr*rmol)
       fac = ((pdz-1.0)/(pzr-1.0))*((pzr+1.0)/(pdz+1.0))
       a = 0.74*dz*alog(fac) + (0.164/rmol)*(pdz-pzr)
    ELSE IF (rmol==0) THEN
       a = 0.74*(dz*alog(dz/zr)-dz+zr)
    ELSE
       a = 0.74*(dz*alog(dz/zr)-dz+zr) + (2.35*rmol)*(dz-zr)**2
    END IF

    fx = a/(vk*ustar*(dz-zr))

    !      CALCULATE THE DEPOSITION VELOCITIY
    DO nss = 1, nspec
       vgtemp(nss) = vgtemp(nss)/(1.0+vgtemp(nss)*fx)
    END DO

    RETURN
  END SUBROUTINE cellvg


  !========================================================================

  SUBROUTINE apply_drydep(m1,m2,m3,ia,iz,ja,jz,V_dep,sclt,M_dep,sclp &
                         ,dens,rtgt,dzt,dt)

    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3
    INTEGER , INTENT(IN)    :: ia
    INTEGER , INTENT(IN)    :: iz
    INTEGER , INTENT(IN)    :: ja
    INTEGER , INTENT(IN)    :: jz

    REAL    , INTENT(IN)    :: V_dep(m2,m3)   ! dry deposition velocity (m/s)
    REAL    , INTENT(INOUT) :: sclt(m1,m2,m3) ! tendency (kg[tracer] kg[air]^-1 s^-1)
    REAL    , INTENT(INOUT) :: M_dep(m2,m3)   ! accumulated mass on surface due dry dep 
                              ! process (kg m^-2)
    REAL    , INTENT(IN)    :: sclp(m1,m2,m3) ! tracer mixing ratio (kg/kg)
    REAL    , INTENT(IN)    :: dens(m2,m3)
    REAL    , INTENT(IN)    :: rtgt(m2,m3) ! (DMK) not used
    REAL    , INTENT(IN)    :: dzt(m1)

    REAL    , INTENT(IN)    :: dt

    INTEGER :: i,j
    REAL    :: dz,tend_drydep			    

    DO j=ja,jz
       DO i=ia,iz
          !-1st vertical layer thickness
          dz = rtgt(i,j)/dzt(2) ! dzt=1/(z(k)-z(k-1))
          !- tendency to dry deposition process	
          tend_drydep  = - V_dep(i,j)*sclp(2,i,j)/(dz)
          !- update total tendency kg[tracer]/kg[air]/s
          sclt(2,i,j) = sclt(2,i,j) + tend_drydep 
          !- accumulate the surface deposited mass of the tracer by this process
          !	M_dep(i,j) = M_dep(i,j) + dt*tend_drydep*dens(i,j) ! kg[tracer] m^-2
          M_dep(i,j) = M_dep(i,j) - dt*tend_drydep*dens(i,j)*dz ! kg[tracer] m^-2
       ENDDO
    ENDDO
    RETURN
  END SUBROUTINE apply_drydep
  !========================================================================

  SUBROUTINE define_PBL_height(m1,m2,m3,npatch,ia,iz,ja,jz,zt,tke,shf,rcp,rtgt,Zi)

    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3
    INTEGER , INTENT(IN)    :: npatch ! (DMK) not used
    INTEGER , INTENT(IN)    :: ia
    INTEGER , INTENT(IN)    :: iz
    INTEGER , INTENT(IN)    :: ja
    INTEGER , INTENT(IN)    :: jz

    REAL    , INTENT(IN)    :: zt(m1)
    REAL    , INTENT(IN)    :: tke(m1,m2,m3)
    REAL    , INTENT(IN)    :: shf(m2,m3)
    REAL    , INTENT(IN)    :: rcp(m1,m2,m3)
    REAL    , INTENT(IN)    :: rtgt(m2,m3)
    REAL    , INTENT(INOUT) :: Zi(m2,m3)

    REAL,PARAMETER :: tkethrsh=0.001       !   tke threshold for PBL height in m2/s2
    !   tkmin    = 5.e-4   minimum TKE in RAMS 
    REAL,PARAMETER :: rcmin=1.e-4          !   min liq water = 0.1 g/kg

    REAL pblht
    INTEGER :: i,j,k


    DO j=ja,jz
       DO i=ia,iz

          Zi(i,j) = 0.
          !- convective layer
          IF(shf(i,j) >= 1.e-8) THEN 

             pblht=0.
             DO k=2,m1-1
                pblht=zt(k)*rtgt(i,j) 
                !                  if(i.ge.10.and.i.le.25.and.j.ge.13.and.j.le.25)
                !     &               print*,'i,j,k,z,pbl=',i,j,k,ztn(k,ngrd),pblht              
                IF( rcp(k,i,j) .GT. rcmin     )GOTO 10 ! dry convective layer  ! (DMK) EXIT
                IF( tke(k,i,j) .LE. tkethrsh  )GOTO 10                         ! (DMK) EXIT
             ENDDO
10           CONTINUE

             Zi(i,j)=pblht
          ENDIF
          !print*,'PBLh',Zi(i,j)

       ENDDO
    ENDDO
    RETURN
  END SUBROUTINE define_PBL_height


  !========================================================================
  SUBROUTINE dep_init(nspecies_chem,dvj)

    INTEGER , INTENT(IN) :: nspecies_chem
    REAL    , INTENT(IN) :: dvj(nspecies_chem)

    ! ..
    ! ..
    ! .. Local Scalars ..
    REAL :: sc
    INTEGER :: iland, iseason, l
    INTEGER :: iprt
    ! ..
    ! .. Local Arrays ..
    REAL :: dat1(nlu,dep_seasons), dat2(nlu,dep_seasons),         &
         dat3(nlu,dep_seasons), dat4(nlu,dep_seasons),         &
         dat5(nlu,dep_seasons), dat6(nlu,dep_seasons),         &
         dat7(nlu,dep_seasons)
    ! ..
    ! .. Data Statements ..
    !     RI for stomatal resistance
    !      data ((ri(ILAND,ISEASON),ILAND=1,nlu),ISEASON=1,dep_seasons)/0.10E+11, &
    DATA ((dat1(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
         0.60E+02, 0.60E+02, 0.60E+02, 0.60E+02, 0.70E+02, 0.12E+03, &
         0.12E+03, 0.12E+03, 0.12E+03, 0.70E+02, 0.13E+03, 0.70E+02, &
         0.13E+03, 0.10E+03, 0.10E+11, 0.80E+02, 0.10E+03, 0.10E+11, &
         0.80E+02, 0.10E+03, 0.10E+03, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.10E+11, 0.10E+11, &
         0.70E+02, 0.25E+03, 0.50E+03, 0.10E+11, 0.10E+11, 0.50E+03, &
         0.10E+11, 0.10E+11, 0.50E+03, 0.50E+03, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.10E+11, &
         0.10E+11, 0.70E+02, 0.25E+03, 0.50E+03, 0.10E+11, 0.10E+11, &
         0.50E+03, 0.10E+11, 0.10E+11, 0.50E+03, 0.50E+03, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.70E+02, 0.40E+03, 0.80E+03, 0.10E+11, &
         0.10E+11, 0.80E+03, 0.10E+11, 0.10E+11, 0.80E+03, 0.80E+03, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.12E+03, 0.12E+03, &
         0.12E+03, 0.12E+03, 0.14E+03, 0.24E+03, 0.24E+03, 0.24E+03, &
         0.12E+03, 0.14E+03, 0.25E+03, 0.70E+02, 0.25E+03, 0.19E+03, &
         0.10E+11, 0.16E+03, 0.19E+03, 0.10E+11, 0.16E+03, 0.19E+03, &
         0.19E+03, 0.10E+11, 0.10E+11, 0.10E+11/
    ! ..
    IF (nlu/=25) THEN
       PRINT *, 'number of land use classifications not correct '
       STOP
    END IF
    IF (dep_seasons/=5) THEN
       PRINT *, 'number of dep_seasons not correct '
       STOP
    END IF

    !     SURFACE RESISTANCE DATA FOR DEPOSITION MODEL OF
    !     M. L. WESELY, ATMOSPHERIC ENVIRONMENT 23 (1989) 1293-1304

    !     Seasonal categories:
    !     1: midsummer with lush vegetation
    !     2: autumn with unharvested cropland
    !     3: late autumn with frost, no snow
    !     4: winter, snow on ground and subfreezing
    !     5: transitional spring with partially green short annuals

    !     Land use types:
    !     USGS type                                Wesely type
    !      1: Urban and built-up land              1
    !      2: Dryland cropland and pasture         2
    !      3: Irrigated cropland and pasture       2
    !      4: Mix. dry/irrg. cropland and pasture  2
    !      5: Cropland/grassland mosaic            2
    !      6: Cropland/woodland mosaic             4
    !      7: Grassland                            3
    !      8: Shrubland                            3
    !      9: Mixed shrubland/grassland            3
    !     10: Savanna                              3, always summer
    !     11: Deciduous broadleaf forest           4
    !     12: Deciduous needleleaf forest          5, autumn and winter modi
    !     13: Evergreen broadleaf forest           4, always summer
    !     14: Evergreen needleleaf forest          5
    !     15: Mixed Forest                         6
    !     16: Water Bodies                         7
    !     17: Herbaceous wetland                   9
    !     18: Wooded wetland                       6
    !     19: Barren or sparsely vegetated         8
    !     20: Herbaceous Tundra                    9
    !     21: Wooded Tundra                        6
    !     22: Mixed Tundra                         6
    !     23: Bare Ground Tundra                   8
    !     24: Snow or Ice                          -, always winter
    !     25: No data                              8


    !     Order of data:
    !      |
    !      |   seasonal category
    !     \|/
    !     ---> landuse USGS type  !kml 
    !     1       2       3       4       5       6      ..... 25    !kml
    !     RLU for outer surfaces in the upper canopy
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
          ri((iland),iseason) = dat1(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          ri(iland,iseason) = dat1(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !      data ((rlu(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
    DATA ((dat2(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
         0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.10E+11, 0.25E+04, 0.20E+04, 0.10E+11, &
         0.25E+04, 0.20E+04, 0.20E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
         0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, 0.90E+04, &
         0.20E+04, 0.40E+04, 0.80E+04, 0.10E+11, 0.90E+04, 0.80E+04, &
         0.10E+11, 0.90E+04, 0.80E+04, 0.80E+04, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
         0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, &
         0.90E+04, 0.20E+04, 0.40E+04, 0.80E+04, 0.10E+11, 0.90E+04, &
         0.80E+04, 0.10E+11, 0.90E+04, 0.80E+04, 0.80E+04, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.20E+04, 0.60E+04, 0.90E+04, 0.10E+11, &
         0.90E+04, 0.90E+04, 0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.40E+04, 0.40E+04, &
         0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, &
         0.20E+04, 0.40E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.30E+04, &
         0.10E+11, 0.40E+04, 0.30E+04, 0.10E+11, 0.40E+04, 0.30E+04, &
         0.30E+04, 0.10E+11, 0.10E+11, 0.10E+11/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rlu((iland),iseason) = dat2(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rlu(iland,iseason) = dat2(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !     RAC for transfer that depends on canopy height and density
    !      data ((rac(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+03, &
    DATA ((dat3(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+04, 0.10E+03, &
         0.10E+03, 0.10E+03, 0.10E+03, 0.20E+04, 0.20E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.00E+00, 0.30E+03, 0.20E+04, 0.00E+00, &
         0.30E+03, 0.20E+04, 0.20E+04, 0.00E+00, 0.00E+00, 0.00E+00, &
         0.10E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+04, &
         0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.15E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.17E+04, 0.00E+00, 0.20E+03, 0.17E+04, &
         0.00E+00, 0.20E+03, 0.17E+04, 0.17E+04, 0.00E+00, 0.00E+00, &
         0.00E+00, 0.10E+03, 0.10E+02, 0.10E+02, 0.10E+02, 0.10E+02, &
         0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+04, &
         0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, 0.00E+00, 0.10E+03, &
         0.15E+04, 0.00E+00, 0.10E+03, 0.15E+04, 0.15E+04, 0.00E+00, &
         0.00E+00, 0.00E+00, 0.10E+03, 0.10E+02, 0.10E+02, 0.10E+02, &
         0.10E+02, 0.10E+04, 0.10E+02, 0.10E+02, 0.10E+02, 0.10E+02, &
         0.10E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, 0.00E+00, &
         0.50E+02, 0.15E+04, 0.00E+00, 0.50E+02, 0.15E+04, 0.15E+04, &
         0.00E+00, 0.00E+00, 0.00E+00, 0.10E+03, 0.50E+02, 0.50E+02, &
         0.50E+02, 0.50E+02, 0.12E+04, 0.80E+02, 0.80E+02, 0.80E+02, &
         0.10E+03, 0.12E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.15E+04, &
         0.00E+00, 0.20E+03, 0.15E+04, 0.00E+00, 0.20E+03, 0.15E+04, &
         0.15E+04, 0.00E+00, 0.00E+00, 0.00E+00/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rac((iland),iseason) = dat3(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rac(iland,iseason) = dat3(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !     RGSS for ground surface  SO2
    !      data ((rgss(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.40E+03, &
    DATA ((dat4(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.40E+03, &
         0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.50E+03, 0.35E+03, &
         0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, 0.50E+03, 0.50E+03, &
         0.50E+03, 0.10E+03, 0.10E+01, 0.10E+01, 0.10E+03, 0.10E+04, &
         0.10E+01, 0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, 0.10E+04, &
         0.40E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.50E+03, &
         0.35E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, 0.50E+03, &
         0.50E+03, 0.50E+03, 0.10E+03, 0.10E+01, 0.10E+01, 0.10E+03, &
         0.10E+04, 0.10E+01, 0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, &
         0.10E+04, 0.40E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, &
         0.50E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.35E+03, 0.50E+03, &
         0.50E+03, 0.50E+03, 0.50E+03, 0.20E+03, 0.10E+01, 0.10E+01, &
         0.20E+03, 0.10E+04, 0.10E+01, 0.20E+03, 0.20E+03, 0.10E+04, &
         0.10E+03, 0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, &
         0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, 0.10E+03, &
         0.10E+03, 0.10E+03, 0.50E+03, 0.10E+03, 0.10E+03, 0.10E+01, &
         0.10E+03, 0.10E+03, 0.10E+04, 0.10E+03, 0.10E+03, 0.10E+03, &
         0.10E+04, 0.10E+03, 0.10E+04, 0.50E+03, 0.15E+03, 0.15E+03, &
         0.15E+03, 0.15E+03, 0.50E+03, 0.35E+03, 0.35E+03, 0.35E+03, &
         0.35E+03, 0.50E+03, 0.50E+03, 0.50E+03, 0.50E+03, 0.20E+03, &
         0.10E+01, 0.10E+01, 0.20E+03, 0.10E+04, 0.10E+01, 0.20E+03, &
         0.20E+03, 0.10E+04, 0.10E+03, 0.10E+04/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rgss((iland),iseason) = dat4(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rgss(iland,iseason) = dat4(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !     RGSO for ground surface  O3
    !      data ((rgso(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.30E+03, &
    DATA ((dat5(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.30E+03, &
         0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.20E+03, 0.20E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
         0.20E+03, 0.30E+03, 0.20E+04, 0.10E+04, 0.30E+03, 0.40E+03, &
         0.10E+04, 0.30E+03, 0.30E+03, 0.40E+03, 0.35E+04, 0.40E+03, &
         0.30E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.20E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
         0.20E+03, 0.20E+03, 0.30E+03, 0.20E+04, 0.80E+03, 0.30E+03, &
         0.40E+03, 0.80E+03, 0.30E+03, 0.30E+03, 0.40E+03, 0.35E+04, &
         0.40E+03, 0.30E+03, 0.15E+03, 0.15E+03, 0.15E+03, 0.15E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.30E+03, 0.20E+04, 0.10E+04, &
         0.30E+03, 0.40E+03, 0.10E+04, 0.30E+03, 0.30E+03, 0.40E+03, &
         0.35E+04, 0.40E+03, 0.60E+03, 0.35E+04, 0.35E+04, 0.35E+04, &
         0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, 0.35E+04, &
         0.35E+04, 0.35E+04, 0.20E+03, 0.35E+04, 0.35E+04, 0.20E+04, &
         0.35E+04, 0.35E+04, 0.40E+03, 0.35E+04, 0.35E+04, 0.35E+04, &
         0.40E+03, 0.35E+04, 0.40E+03, 0.30E+03, 0.15E+03, 0.15E+03, &
         0.15E+03, 0.15E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, &
         0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.20E+03, 0.30E+03, &
         0.20E+04, 0.10E+04, 0.30E+03, 0.40E+03, 0.10E+04, 0.30E+03, &
         0.30E+03, 0.40E+03, 0.35E+04, 0.40E+03/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rgso((iland),iseason) = dat5(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rgso(iland,iseason) = dat5(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !     RCLS for exposed surfaces in the lower canopy  SO2
    !      data ((rcls(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
    DATA ((dat6(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
         0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.20E+04, &
         0.20E+04, 0.20E+04, 0.10E+11, 0.25E+04, 0.20E+04, 0.10E+11, &
         0.25E+04, 0.20E+04, 0.20E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, &
         0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, 0.90E+04, &
         0.20E+04, 0.20E+04, 0.40E+04, 0.10E+11, 0.90E+04, 0.40E+04, &
         0.10E+11, 0.90E+04, 0.40E+04, 0.40E+04, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.90E+04, 0.90E+04, 0.90E+04, 0.90E+04, 0.20E+04, 0.90E+04, &
         0.90E+04, 0.20E+04, 0.30E+04, 0.60E+04, 0.10E+11, 0.90E+04, &
         0.60E+04, 0.10E+11, 0.90E+04, 0.60E+04, 0.60E+04, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.90E+04, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.90E+04, 0.90E+04, 0.20E+04, 0.20E+03, 0.40E+03, 0.10E+11, &
         0.90E+04, 0.40E+03, 0.10E+11, 0.90E+04, 0.40E+03, 0.40E+03, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.40E+04, 0.40E+04, &
         0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, 0.40E+04, &
         0.20E+04, 0.40E+04, 0.20E+04, 0.20E+04, 0.20E+04, 0.30E+04, &
         0.10E+11, 0.40E+04, 0.30E+04, 0.10E+11, 0.40E+04, 0.30E+04, &
         0.30E+04, 0.10E+11, 0.10E+11, 0.10E+11/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rcls((iland),iseason) = dat6(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rcls(iland,iseason) = dat6(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO
    !     RCLO for exposed surfaces in the lower canopy  O3
    !      data ((rclo(ILAND,ISEASON),ILAND=1,25),ISEASON=1,5)/0.10E+11, &
    DATA ((dat7(iland,iseason),iland=1,nlu),iseason=1,dep_seasons)/0.10E+11, &
         0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
         0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
         0.10E+04, 0.10E+04, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+11, &
         0.10E+04, 0.10E+04, 0.10E+04, 0.10E+11, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, &
         0.40E+03, 0.40E+03, 0.40E+03, 0.10E+04, 0.40E+03, 0.40E+03, &
         0.10E+04, 0.10E+04, 0.60E+03, 0.10E+11, 0.40E+03, 0.60E+03, &
         0.10E+11, 0.40E+03, 0.60E+03, 0.60E+03, 0.10E+11, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
         0.40E+03, 0.40E+03, 0.40E+03, 0.40E+03, 0.10E+04, 0.40E+03, &
         0.40E+03, 0.10E+04, 0.10E+04, 0.60E+03, 0.10E+11, 0.80E+03, &
         0.60E+03, 0.10E+11, 0.80E+03, 0.60E+03, 0.60E+03, 0.10E+11, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, 0.10E+04, &
         0.10E+04, 0.40E+03, 0.10E+04, 0.10E+04, 0.10E+04, 0.10E+04, &
         0.40E+03, 0.40E+03, 0.10E+04, 0.15E+04, 0.60E+03, 0.10E+11, &
         0.80E+03, 0.60E+03, 0.10E+11, 0.80E+03, 0.60E+03, 0.60E+03, &
         0.10E+11, 0.10E+11, 0.10E+11, 0.10E+11, 0.10E+04, 0.10E+04, &
         0.10E+04, 0.10E+04, 0.50E+03, 0.50E+03, 0.50E+03, 0.50E+03, &
         0.10E+04, 0.50E+03, 0.15E+04, 0.10E+04, 0.15E+04, 0.70E+03, &
         0.10E+11, 0.60E+03, 0.70E+03, 0.10E+11, 0.60E+03, 0.70E+03, &
         0.70E+03, 0.10E+11, 0.10E+11, 0.10E+11/
    DO iseason = 1, dep_seasons
       DO iland = 1, nlu

!--(DMK-CCATT-INI)----------------------------------------------------------
         rclo((iland),iseason) = dat7(iland,iseason)
!--(DMK-CCATT-OLD)----------------------------------------------------------
!          rclo(iland,iseason) = dat7(iland,iseason)
!--(DMK-CCATT-FIM)----------------------------------------------------------

       END DO
    END DO

    DO l = 1, nspecies_chem
       ! if(dvj(l) < 0.) dvj(l) = 1./(sqrt(weight(l)))
       !          hstar4(l) = hstar(l) ! preliminary              
       ! Correction of diff. coeff
       !          dvj(l) = dvj(l)*(293.15/298.15)**1.75
       !          sc = 0.15/dvj(l) ! Schmidt Number at 20�C   
       dratio(l) = 0.242/dvj(l) !    ! of water vapor and gas at
       ! Ratio of diffusion coeffi
       !          scpr23(l) = (sc/0.72)**(2./3.) ! (Schmidt # / Prandtl #)**
    END DO


    !     DATA FOR AEROSOL PARTICLE DEPOSITION FOR THE MODEL OF
    !     J. W. ERISMAN, A. VAN PUL AND P. WYERS
    !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607

    !     vd = (u* / k) * CORRECTION FACTORS

    !     CONSTANT K FOR LANDUSE TYPES:
    !srf- general initialization
    kpart(:) = 500.
    ! urban and built-up land                  
    kpart(19) = 500.
    ! dryland cropland and pasture             
    kpart(2) = 500.
    ! irrigated cropland and pasture           
    kpart(16) = 500.
    ! mixed dryland/irrigated cropland and past
    kpart(16) = 500.
    ! cropland/grassland mosaic                
    kpart(15) = 500.
    ! cropland/woodland mosaic                 
    kpart(14) = 100.
    ! grassland                                
    kpart(8) = 500.
    ! shrubland                                
    kpart(13) = 500.
    ! mixed shrubland/grassland                
    kpart(9) = 500.
    ! savanna                                  
    kpart(9) = 500.
    ! deciduous broadleaf forest               
    kpart(5) = 100.
    ! deciduous needleleaf forest              
    kpart(6) = 100.
    ! evergreen broadleaf forest               
    kpart(7) = 100.
    ! evergreen needleleaf forest              
    kpart(4) = 100.
    ! mixed forest                             
    kpart(14) = 100.
    ! water bodies                             
    kpart(1) = 500.
    ! herbaceous wetland                       
    kpart(17) = 500.
    ! wooded wetland                           
    kpart(18) = 500.
    ! barren or sparsely vegetated             
    kpart(10) = 500.
    ! herbaceous tundra                        
    kpart(11) = 500.
    ! wooded tundra                            
    kpart(11) = 100.
    ! mixed tundra                             
    kpart(11) = 500.
    ! bare ground tundra                       
    kpart(3) = 500.
    ! snow or ice                              
    kpart(2) = 500.
    !     Comments:
    kpart(25) = 500.
    !     Erisman et al. (1994) give
    !     k = 500 for low vegetation and k = 100 for forests.

    !     For desert k = 500 is taken according to measurements
    !     on bare soil by
    !     J. Fontan, A. Lopez, E. Lamaud and A. Druilhet (1997)
    !     "Vertical Flux Measurements of the Submicronic Aerosol Particles
    !     and Parametrisation of the Dry Deposition Velocity"
    !     in: "Biosphere-Atmosphere Exchange of Pollutants
    !     and Trace Substances"
    !     Editor: S. Slanina. Springer-Verlag Berlin, Heidelberg, 1997
    !     pp. 381-390

    !     For coniferous forest the Erisman value of  k = 100 is taken.
    !     Measurements of Erisman et al. (1997) in a coniferous forest
    !     in the Netherlands, lead to values of k between 20 and 38
    !     (Atmospheric Environment 31 (1997), 321-332).
    !     However, these high values of vd may be reached during
    !     instable cases. The eddy correlation measurements
    !     of Gallagher et al. (1997) made during the same experiment
    !     show for stable cases (L>0) values of k between 200 and 250
    !     at minimum (Atmospheric Environment 31 (1997), 359-373).
    !     Fontan et al. (1997) found k = 250 in a forest
    !     of maritime pine in southwestern France.

    !     For gras, model calculations of Davidson et al. support
    !     the value of 500.
    !     C. I. Davidson, J. M. Miller and M. A. Pleskov
    !     "The Influence of Surface Structure on Predicted Particles
    !     Dry Deposition to Natural Gras Canopies"
    !     Water, Air, and Soil Pollution 18 (1982) 25-43

    !     Snow covered surface: The experiment of Ibrahim et al. (1983)
    !     gives k = 436 for 0.7 um diameter particles.
    !     The deposition velocity of Milford and Davidson (1987)
    !     gives k = 154 for continental sulfate aerosol.
    !     M. Ibrahim, L. A. Barrie and F. Fanaki
    !     Atmospheric Environment 17 (1983), 781-788

    !     J. B. Milford and C. I. Davidson
    !     "The Sizes of Particulate Sulfate and Nitrate in the Atmosphere
    !     - A Review"
    !     JAPCA 37 (1987), 125-134
    !
    !------------!-------------!----------------------------------------------------------------
    !  LEAF-3 CLASS (20)                      !  Wesely type !   Land use types:
    !---AND DESCRIPTION 			  !		 !   USGS type  			 
    !  0  Ocean                               !   7 	 !    1: Urban and built-up land	      
    !  1  Lakes, rivers, streams              !   7 	 !    2: Dryland cropland and pasture	     
    !  2  Ice cap/glacier                     !   - 	 !    3: Irrigated cropland and pasture      
    !  3  Desert, bare soil                   !   8 	 !    4: Mix. dry/irrg. cropland and pasture 
    !  4  Evergreen needleleaf tree           !   5 	 !    5: Cropland/grassland mosaic	     
    !  5  Deciduous needleleaf tree           !   5 	 !    6: Cropland/woodland mosaic	     
    !  6  Deciduous broadleaf tree            !   4 	 !    7: Grassland			     
    !  7  Evergreen broadleaf tree            !   4 	 !    8: Shrubland			     
    !  8  Short grass                         !   2 	 !    9: Mixed shrubland/grassland				
    !  9  Tall grass                          !   3 	 !   10: Savanna			     
    ! 10  Semi-desert                         !   8 	 !   11: Deciduous broadleaf forest	     
    ! 11  Tundra                              !   8 	 !   12: Deciduous needleleaf forest	     
    ! 12  Evergreen shrub                     !   3 	 !   13: Evergreen broadleaf forest	     
    ! 13  Deciduous shrub                     !   3 	 !   14: Evergreen needleleaf forest	     
    ! 14  Mixed woodland                      !   6 	 !   15: Mixed Forest			     
    ! 15  Crop/mixed farming, C3 grassland    !   2 	 !   16: Water Bodies			     
    ! 16  Irrigated crop                      !   2 	 !   17: Herbaceous wetland		     
    ! 17  Bog or marsh                        !   9 	 !   18: Wooded wetland 		     
    ! 18  Wooded grassland                    !   6 	 !   19: Barren or sparsely vegetated	     
    ! 19  Urban and built up                  !   1 	 !   20: Herbaceous Tundra		     
    ! 20  Wetland evergreen broadleaf tree    !   6 	 !   21: Wooded Tundra  		     
    !		                                         !   22: Mixed Tundra			     
    !		                                         !   23: Bare Ground Tundra		     
    !	 	                                         !   24: Snow or Ice			       
    !--------------------------------------------------------------------------------------------
!KML: Os valores de ixxxlu foram trocados para fazer a equivalencia entre as classes do
! do LEAF3 e do USGS porque as tabelas do Wesely estao em funcao das classes do USGS (e nao das classes do
!Wesely).

!kml    ixxxlu(1)  =7 
!kml    ixxxlu(2)  =7 
!kml    ixxxlu(3)  =8 
!kml    ixxxlu(4)  =5 
!kml    ixxxlu(5)  =5 
!kml    ixxxlu(6)  =4 
!kml    ixxxlu(7)  =4 
!kml    ixxxlu(8)  =2 
!kml    ixxxlu(9)  =3 
!kml    ixxxlu(10) =8 
!kml    ixxxlu(11) =8 
!kml    ixxxlu(12) =3 
!kml    ixxxlu(13) =3 
!kml    ixxxlu(14) =6 
!kml    ixxxlu(15) =2 
!kml    ixxxlu(16) =2 
!kml    ixxxlu(17) =9 
!kml    ixxxlu(18) =6 
!kml    ixxxlu(19) =1 
!kml    ixxxlu(20) =6 



    ixxxlu(1)  =16 
    ixxxlu(2)  =24 
    ixxxlu(3)  =19 
    ixxxlu(4)  =14 
    ixxxlu(5)  =12 
    ixxxlu(6)  =11 
    ixxxlu(7)  =13 
    ixxxlu(8)  =7 
    ixxxlu(9)  =8 
    ixxxlu(10) =19 
    ixxxlu(11) =22 
    ixxxlu(12) =8 
    ixxxlu(13) =8 
    ixxxlu(14) =15 
    ixxxlu(15) =5 
    ixxxlu(16) =3 
    ixxxlu(17) =17 
    ixxxlu(18) =10 
    ixxxlu(19) =1 
    ixxxlu(20) =18 


  END SUBROUTINE dep_init

  ! **********************************************************************
  SUBROUTINE deppart(rmol,ustar,rh,clw,iland,dvpart,dvfog)
    !     THIS SUBROUTINE CALCULATES SURFACE DEPOSITION VELOCITIES
    !     FOR FINE AEROSOL PARTICLES ACCORDING TO THE MODEL OF
    !     J. W. ERISMAN, A. VAN PUL, AND P. WYERS,
    !     ATMOSPHERIC ENVIRONMENT 28 (1994), 2595-2607
    !     WRITTEN BY WINFRIED SEIDL, APRIL 1997
    !     MODIFIED BY WINFRIED SEIDL, MARCH 2000
    !            FOR MM5 VERSION 3
    ! ----------------------------------------------------------------------

    REAL    , INTENT(IN)  :: rmol
    REAL    , INTENT(IN)  :: ustar
    REAL    , INTENT(IN)  :: rh
    REAL    , INTENT(IN)  :: clw
    INTEGER , INTENT(IN)  :: iland
    REAL    , INTENT(OUT) :: dvpart
    REAL    , INTENT(OUT) :: dvfog

    ! ..
    ! .. Intrinsic Functions ..
    INTRINSIC exp
    ! ..
    dvpart = ustar/kpart(iland)
    !print*,ustar,iland,kpart(iland)

    IF (rmol<0.) THEN
       !         INSTABLE LAYERING CORRECTION
       dvpart = dvpart*(1.+(-300.*rmol)**0.66667)
    END IF
    IF (rh>80.) THEN
       !         HIGH RELATIVE HUMIDITY CORRECTION
       !         ACCORDING TO J. W. ERISMAN ET AL.
       !         ATMOSPHERIC ENVIRONMENT 31 (1997), 321-332
       dvpart = dvpart*(1.+0.37*EXP((rh-80.)/20.))
    END IF
    !print*,'dvpart=', dvpart
    !srf - not using fog deposition (in case on, set clw,iland=iveg)
    RETURN

    !       SEDIMENTATION VELOCITY OF FOG WATER ACCORDING TO
    !       R. FORKEL, W. SEIDL, R. DLUGI AND E. DEIGELE
    !       J. GEOPHYS. RES. 95D (1990), 18501-18515
    dvfog = 0.06*clw
!kml    IF (ixxxlu(iland)==5) THEN
    IF (ixxxlu(iland)==12 .OR.ixxxlu(iland)==14 ) THEN
       !         TURBULENT DEPOSITION OF FOG WATER IN CONIFEROUS FOREST ACCORDI
       !         A. T. VERMEULEN ET AL.
       !         ATMOSPHERIC ENVIRONMENT 31 (1997), 375-386
       dvfog = dvfog + 0.195*ustar*ustar
    END IF

  END SUBROUTINE deppart
  ! **********************************************************************
  SUBROUTINE depvel(numchem,rmol,zr,z0,ustar,vgpart)!,polint)
    !     THIS FUNCTION HAS BEEN DESIGNED TO EVALUATE AN UPPER LIMIT
    !     FOR THE POLLUTANT DEPOSITION VELOCITY AS A FUNCTION OF THE
    !     SURFACE ROUGHNESS AND METEOROLOGICAL CONDITIONS.
    !     PROGRAM WRITTEN BY GREGORY J.MCRAE (NOVEMBER 1977)
    !         Modified by Darrell A. Winner  (Feb. 1991)
    !                  by Winfried Seidl     (Aug. 1997)
    !.....PROGRAM VARIABLES...
    !     RMOL     - RECIPROCAL OF THE MONIN-OBUKHOV LENGTH
    !     ZR       - REFERENCE HEIGHT
    !     Z0       - SURFACE ROUGHNESS HEIGHT
    !     SCPR23   - (Schmidt #/Prandtl #)**(2/3) Diffusion correction fact
    !     UBAR     - ABSOLUTE VALUE OF SURFACE WIND SPEED
    !     DEPVEL   - POLLUTANT DEPOSITION VELOCITY
    !     Vk       - VON KARMAN CONSTANT
    !     USTAR    - FRICTION VELOCITY U*
    !     POLINT   - POLLUTANT INTEGRAL
    !.....REFERENCES...
    !     MCRAE, G.J. ET AL. (1983) MATHEMATICAL MODELING OF PHOTOCHEMICAL
    !       AIR POLLUTION, ENVIRONMENTAL QUALITY LABORATORY REPORT 18,
    !       CALIFORNIA INSTITUTE OF TECHNOLOGY, PASADENA, CALIFORNIA.
    !.....RESTRICTIONS...
    !     1. THE MODEL EDDY DIFFUSIVITIES ARE BASED ON MONIN-OBUKHOV
    !        SIMILARITY THEORY AND SO ARE ONLY APPLICABLE IN THE
    !        SURFACE LAYER, A HEIGHT OF O(30M).
    !     2. ALL INPUT UNITS MUST BE CONSISTENT
    !     3. THE PHI FUNCTIONS USED TO CALCULATE THE FRICTION
    !        VELOCITY U* AND THE POLLUTANT INTEGRALS ARE BASED
    !        ON THE WORK OF BUSINGER ET AL.(1971).
    !     4. THE MOMENTUM AND POLLUTANT DIFFUSIVITIES ARE NOT
    !        THE SAME FOR THE CASES L<0 AND L>0.

    INTEGER , INTENT(IN)  :: numchem ! (DMK) not used
    REAL    , INTENT(IN)  :: rmol
    REAL    , INTENT(IN)  :: zr
    REAL    , INTENT(IN)  :: z0
    REAL    , INTENT(IN)  :: ustar
    REAL    , INTENT(OUT) :: vgpart

    ! ..
    ! .. Array Arguments ..
    !        REAL :: depv(numchem)
    ! ..
    ! .. Local Scalars ..
    REAL :: ao, ar, polint
    INTEGER :: l
    ! ..
    ! .. Intrinsic Functions ..
    INTRINSIC alog
    ! ..
    !     Set the von Karman constant
    REAL, PARAMETER :: vk = 0.4

    !     DETERMINE THE STABILITY BASED ON THE CONDITIONS
    !             1/L < 0 UNSTABLE
    !             1/L = 0 NEUTRAL
    !             1/L > 0 STABLE

    IF (rmol<0) THEN
       ar = ((1.0-9.0*zr*rmol)**(0.25)+0.001)**2
       ao = ((1.0-9.0*z0*rmol)**(0.25)+0.001)**2
       polint = 0.74*(alog((ar-1.0)/(ar+1.0))-alog((ao-1.0)/(ao+1.0)))
    ELSE IF (rmol==0) THEN
       polint = 0.74*alog(zr/z0)
    ELSE
       polint = 0.74*alog(zr/z0) + 4.7*rmol*(zr-z0)
    END IF

    !     CALCULATE THE Maximum DEPOSITION VELOCITY

    !DO l = 1, numchem
    !  depv(l) = ustar*vk/(2.0*scpr23(l)+polint)
    !END DO
    vgpart = ustar*vk/polint
    !print*,polint,ustar,vgpart,z0,ustar*vk/polint
    RETURN
  END SUBROUTINE depvel

  !========================================================================

  ! -----------------   AEROSOL DRY DEP AND SEDIMENTATION -----------------

  !========================================================================

  SUBROUTINE dry_dep_sedim_particles(m1,m2,m3,naddsc,npatch              &
                                    ,ia,iz,ja,jz,r_aer,temp3d,air_dens3d &
                                    ,temps,dens,vels,rvs,Zi,ustar,tstar  &
                                    ,patch_area,veg,Z0m,nmodes,nspecies  &
                                    ,part_radius,part_dens,mode_alloc,on &
                                    ,g,dd_sedim)

    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: naddsc ! (DMK) not used
    INTEGER , INTENT(IN) :: npatch
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz

    REAL    , INTENT(IN) :: r_aer(m2,m3,npatch)
    REAL    , INTENT(IN) :: temp3d(m1,m2,m3)
    REAL    , INTENT(IN) :: air_dens3d(m1,m2,m3)
    REAL    , INTENT(IN) :: temps(m2,m3)
    REAL    , INTENT(IN) :: dens(m2,m3)
    REAL    , INTENT(IN) :: vels(m2,m3)
    REAL    , INTENT(IN) :: rvs(m2,m3) ! (DMK) not used
    REAL    , INTENT(IN) :: Zi(m2,m3)
    REAL    , INTENT(IN) :: ustar(m2,m3,npatch)
    REAL    , INTENT(IN) :: tstar(m2,m3,npatch)
    REAL    , INTENT(IN) :: patch_area(m2,m3,npatch)
    REAL    , INTENT(IN) :: veg(m2,m3,npatch)
    REAL    , INTENT(IN) :: Z0m(m2,m3,npatch)

    ! aer1_list
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: nspecies
    REAL    , INTENT(IN) :: part_radius(nmodes,nspecies)
    REAL    , INTENT(IN) :: part_dens(nmodes,nspecies)
    INTEGER , INTENT(IN) :: mode_alloc(nmodes,nspecies_aer)
    INTEGER , INTENT(IN) :: on

    ! rconstants
    REAL    , INTENT(IN) :: g

    TYPE(sedim_type), INTENT(INOUT) :: dd_sedim(naer_transported)

    REAL    :: vdtmp
    INTEGER :: i,j,ipatch,ispc

!    IF(.NOT. aer_alloc) THEN
!       CALL alloc_aer_sedim(m1,m2,m3,npatch,ngrids, &
!                            nmodes,nspecies,mode_alloc, &
!                            on,mmzp,mmxp,mmyp) 
!    END IF

    IF( NAER_TRANSPORTED == 0 ) RETURN

    !- sedimentation  parameterization 
    CALL sedim_particles_3d(m1,m2,m3,npatch,ia,iz,ja,jz,temp3d,air_dens3d, &
                            nmodes,nspecies,part_radius,part_dens,g,dd_sedim)


    !- dry deposition parameterization

    !- laminar sub-layer resistance
    CALL lsl_particles(m2,m3,npatch,ia,iz,ja,jz &
         ,temps,dens,vels,rvs,Zi,ustar,tstar,patch_area,veg,Z0m &
         ,nmodes,nspecies,part_radius,part_dens &
         ,dd_sedim)
    
    !- particles deposition velocity (m/s)

    DO j = ja,jz
       DO i = ia,iz
          DO ispc=naer_a,naer_z

             dd_sedim(ispc)%v_dep_part(i,j) = 0.

             DO ipatch = 1,npatch

                IF (patch_area(i,j,ipatch) .GE. .009) THEN


                   vdtmp = dd_sedim(ispc)%v_sed_part(1,i,j) + &
                        1./( r_aer(i,j,ipatch) + dd_sedim(ispc)%r_lsl_part(i,j,ipatch) + &
                        r_aer(i,j,ipatch) * dd_sedim(ispc)%r_lsl_part(i,j,ipatch) * &
                        dd_sedim(ispc)%v_sed_part(1,i,j)				  ) 

                   dd_sedim(ispc)%v_dep_part(i,j) = dd_sedim(ispc)%v_dep_part(i,j)  + &
                        patch_area(i,j,ipatch)*vdtmp

                ENDIF
             ENDDO

             !print*,'V_dep (cm/s)',i,j,dd_sedim(ispc,ngrid)%v_dep_part(i,j)*100.
             !print*,'r_lsl-v_sed =', dd_sedim(ispc,ngrid)%r_lsl_part(i,j,1:npatch), &
             !	dd_sedim(ispc,ngrid)%v_sed_part(1,i,j) !,r_aer(i,j,ipatch)

          ENDDO
       ENDDO
    ENDDO

    RETURN

    !- replace the sedimentation velocity at surface with the deposition velocity
    !- (which includes v_sed_part plus the turbulent flux in laminar sub-layer)

    DO j = ja,jz
       DO i = ia,iz
          DO ispc=naer_a,naer_z
             dd_sedim(ispc)%v_sed_part(1,i,j) = dd_sedim(ispc)%v_dep_part(i,j)
          ENDDO
       ENDDO
    ENDDO


  END SUBROUTINE dry_dep_sedim_particles
  !========================================================================

  SUBROUTINE lsl_particles(m2,m3,npatch,ia,iz,ja,jz  &
                          ,temps,dens,vels,rvs,Zi,ustar &
                          ,tstar,patch_area,veg,Z0m     &
                          ,nmodes,nspecies,part_radius  &
                          ,part_dens                    &
                          ,dd_sedim)

    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: npatch
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz

    REAL    , INTENT(IN) :: temps(m2,m3)
    REAL    , INTENT(IN) :: dens(m2,m3)
    REAL    , INTENT(IN) :: vels(m2,m3)
    REAL    , INTENT(IN) :: rvs(m2,m3) ! (DMK) not used
    REAL    , INTENT(IN) :: Zi(m2,m3)
    REAL    , INTENT(IN) :: ustar(m2,m3,npatch)
    REAL    , INTENT(IN) :: tstar(m2,m3,npatch)
    REAL    , INTENT(IN) :: patch_area(m2,m3,npatch)
    REAL    , INTENT(IN) :: veg(m2,m3,npatch)
    REAL    , INTENT(IN) :: Z0m(m2,m3,npatch) ! (DMK) not used

    ! aer1_list
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: nspecies
    REAL    , INTENT(IN) :: part_radius(nmodes,nspecies)
    REAL    , INTENT(IN) :: part_dens(nmodes,nspecies)

    TYPE(sedim_type), INTENT(INOUT) :: dd_sedim(naer_transported)

    REAL,PARAMETER :: kB = 1.3807e-23      ! const Boltzmann - kg m^2 s^-2 K^-1 molecule^-1
    REAL,PARAMETER :: ASP = 1.257          ! 1.249
    REAL,PARAMETER :: BSP = 0.4            ! 0.42
    REAL,PARAMETER :: CSP = 1.1            ! 0.87 

    REAL,PARAMETER :: M_AVEG = 4.8096e-26  ! average mass of one molecure - kg  molecule^-1
    REAL,PARAMETER :: vonK = 0.40          ! von Karman constant
    REAL,PARAMETER :: Cpd = 1004.67        ! specific heat of dry air [J/kg/K]
    REAL,PARAMETER :: em23 = -2./3., ep23 = +2./3., ep13=1./3.    ! exponents 2./3. 1/3
    REAL,PARAMETER :: pi = 3.1415927, g = 9.80

    INTEGER :: i,j,ipatch,ispc
    REAL    :: wptp,wstar
    !
    ! dd_sedim(ispc,ng)%v_sed      == particle sedimentation velocity   (m/s)
    ! dd_sedim(ispc,ng)%r_lsl_part == resistance to molecular diffusion (s/m)
    !
    REAL Kn,n_air,Gi,v_air  ,mfp,D,nu,Sc,St,Kd,Dh,Z0h,Pr    

    REAL,PARAMETER :: limite = -30. 

    DO j = ja,jz
       DO i = ia,iz

          !- several particle/environment properties

          !- mean speed of air molecules (m/s)
          !  v_air = sqrt( 8. * kB   * temps(i,j) / (pi * M_AVEG) )
          v_air = SQRT( 7.3102e+2 * temps(i,j)			  )

          !-dynamic viscosity of air (kg m^-1 s^-1)
          !  n_air = 1.8325e-5*(416.16/(temps(i,j)+120.))*(temps(i,j)/296.16)**1.5
          !optimized version
          n_air = 1.8325e-5*(416.16/(temps(i,j)+120.))*(temps(i,j)/296.16) * & 
               SQRT(temps(i,j)/296.16)

          !- mean free path of an air molecule (m)
          mfp = 2.* n_air /(dens(i,j)*v_air)

          !-  kinematic viscosity of air ()
          nu = n_air/dens(i,j)

          DO ispc=naer_a,naer_z

             !- Knudsen number
             Kn = mfp/part_radius(ind_mode(ispc),ind_aer(ispc))


             !- Slip correction factor (Gi)
             Gi = 1. + Kn*( ASP + BSP*EXP(-CSP/Kn) )

             !- Schmidt number determination (Jacobson)
             !-- molecular diffusivity (Brownian diffusivity coeficient)   
             !    D =  (kB/M_AVEG) * temps(i,j) * Gi / (6.*pi*part_radius*n_air)

             D =  kB * temps(i,j) * Gi / (6.*pi*part_radius(ind_mode(ispc),ind_aer(ispc))*n_air)


             !-  Schmidt number
             Sc = nu/D

             !- laminar sub-layer resistance for particles (s m^-1)
             DO ipatch=1,npatch

                !-  Stokes number determination (Slinn, 1980)
                !-  St =     ustar^2 * V_sedim /( g * kinematic viscosity of air)    
                St = MAX(ustar(i,j,ipatch)**2. * dd_sedim(ispc)%v_sed_part(1,i,j) &
                     / (g * nu), 0.01) 

                !-  laminar sub-layer resistance for particles (s m^-1) 

                IF(ipatch == 1) THEN !- water

                   !-from Slinn 1980 (Atmos. Env.)
                   !dd_sedim(ispc,ng)%r_lsl_part(i,j,ipatch) = ( vonK * vels(i,j)/ustar(i,j,ipatch)**2. ) / &
                   !	                  ( 1./sqrt(Sc) + 10.**(-3./St) )
                   dd_sedim(ispc)%r_lsl_part(i,j,ipatch) = ( vonK * vels(i,j)/ustar(i,j,ipatch)**2. ) / &
                        ( 1./SQRT(Sc) + 10.**MAX((-3./St),limite) )


                ELSE
                   IF(patch_area(i,j,ipatch) > 0.009) THEN

                      !-- for smooth land surfaces !- bare ground (3), desert(3) or ice (2)
                      !-- Seinfeld & Pandis (1998)
                      IF(NINT(veg(i,j,ipatch)) == 3 ) THEN 

                         ! dd_sedim(ispc,ng)%r_lsl_part(i,j,ipatch) = 1./(ustar(i,j,ipatch) * (Sc**em23 + 10.**(-3./St)))
                         dd_sedim(ispc)%r_lsl_part(i,j,ipatch) = 1./(ustar(i,j,ipatch) * &
                              (Sc**em23 + 10.**MAX((-3./St),limite)))


                      ELSE
                         !- lsl resistance according  Jacobson(1999)
                         !
                         !- thermal conductivity of dry air (Kd)
                         !Kd = 0.023807 + 7.1128e-5*(temps(i,j) - 273.15) !- Eq.(2.3) 
                         !- Prandt number
                         !Pr =  n_air*Cpd*(1.+0.859*rvs(i,j))/Kd           !- Eq.(17.32)  

                         !- energy moisture roughness lengths (Z0h)                 !- Eq.(8.10)
                         !-- molecular thermal diffusion coeff. (m^2 s^-1)
                         !Dh=Kd/(dens(i,j)*Cpd)
                         !- Z0h	   
                         !Z0h=Dh/(vonK*ustar(i,j,ipatch))
                         !- lsl resistance according Jacobson Eq. (20.14)
                         ! dd_sedim(ispc,ng)%r_lsl_part(i,j,ipatch) =  log( Z0m(i,j,ipatch)/Z0h )* ( (Sc/Pr)**ep23 ) &
                         !	                     /  ( vonK*ustar(i,j,ipatch) )

                         !---------  		 
                         !- lsl resistance according :
                         !- from Wesely et al. (1985) following Slinn (1982)
                         !- also Binkowski & Shankar, JGR, 100,D12,26191-26209, 1995
                         !- Rb= (u* (1+ (w*/u*)^2)) (Sc^2./3 + 10^(-3/St))	 

                         wptp  = - ustar(i,j,ipatch) * tstar(i,j,ipatch) ! sensible heat flux
                         wstar = ( MAX (0., g* Zi(i,j)*  wptp/temps(i,j) ) )**ep13  
                         !	     
                         !dd_sedim(ispc,ng)%r_lsl_part(i,j,ipatch) = 1./(  					  &
                         !				     ustar(i,j,ipatch)* 			 &
                         !				     (1. + 0.24*(wstar/ustar(i,j,ipatch))**2.)*  &
                         !				     ( (Sc**em23 + 10.**(-3./St)) )		 &
                         !				    )
                         dd_sedim(ispc)%r_lsl_part(i,j,ipatch) = 1./(						&
                              ustar(i,j,ipatch)*  			&
                              (1. + 0.24*(wstar/ustar(i,j,ipatch))**2.)*  &
                              ( (Sc**em23 + 10.**MAX((-3./St),limite) ) ) )

                         !       print*,'================'
                         !       print*,'LSL_WE',i,j,ipatch,nint(veg(i,j,ipatch)),patch_area(i,j,ipatch)
                         !	PRINT*,wstar,ustar(i,j,ipatch),r_lsl(i,j,ipatch)

                      ENDIF
                   ENDIF
                ENDIF


                !print*,'laminar sub-layer resistance'
                !print*,'LSL_J',i,j,ipatch,dd_sedim(ispc,ng)%r_lsl_part(i,j,ipatch)
                !print*,'ZOM',Z0m(i,j,ipatch)
                !print*,'Z0H',Z0h
                !print*,'SC PR',Sc,Pr
                !print*,'u*',ustar(i,j,ipatch)

             ENDDO
          ENDDO

       ENDDO
    ENDDO

  END SUBROUTINE lsl_particles


  !------------------------------------------------------------------------
  SUBROUTINE sedim_particles_3d(m1,m2,m3,npatch,ia,iz,ja,jz,temp3d       &
                               ,air_dens3d,nmodes,nspecies,part_radius   &
                               ,part_dens,g                              &
                               ,dd_sedim)

    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: npatch
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    REAL    , INTENT(IN) :: temp3d(m1,m2,m3)
    REAL    , INTENT(IN) :: air_dens3d(m1,m2,m3)

    ! aer1_list
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: nspecies
    REAL    , INTENT(IN) :: part_radius(nmodes,nspecies)
    REAL    , INTENT(IN) :: part_dens(nmodes,nspecies)

    ! rconstants
    REAL    , INTENT(IN) :: g

    TYPE(sedim_type), INTENT(INOUT) :: dd_sedim(naer_transported)

    REAL,PARAMETER :: ASP = 1.257          ! 1.249
    REAL,PARAMETER :: BSP = 0.4            ! 0.42
    REAL,PARAMETER :: CSP = 1.1            ! 0.87 
    REAL,PARAMETER :: onesix = 1./6.         

    INTEGER :: i,j,ipatch,k,ispc

    !
    !For small particles diameter (<20 micrometer) (low Reynolds number)
    !the fall velocity is given by (Seinfeld & Pandis, Jacobson)
    ! V_s = 2 r**2 (rho_p - rho_air) * g * Gi / 9 n_air
    ! where:
    ! r == radius of particle (m)
    ! rho_p,a = density of particle, air (kg/m^3)
    ! g = gravity accel (9.8 m/s^2)
    ! Gi = Kn(A' + B' + C' exp(-C'/Kn) , Kn = Knudsen number
    ! n_air = dynamic viscosity of air
    !- constantes
    !parameter (ASP=1.257,BSP=0.4,CSP=1.1) ! Seinfeld & Pandis 
    !parameter (kB = 1.3807e-23)     ! const Boltzmann - kg m^2 s^-2 K^-1 molecule^-1
    !parameter (M_AVEG = 4.8096e-26) ! average mass of one molecure - kg  molecule^-1
    !
    !real, dimension(m1,m2,m3,nmodes,aer_nspecies)    :: v_sed

    REAL A,B,C,Kn,n_air,Gi,v_air,mfp,nu,re,z,b0,x,y
    !- for FALL 3d formulation
    REAL Cd
    REAL, PARAMETER :: eps = 1.e-6
    REAL, PARAMETER :: psi = 0.95
    REAL, PARAMETER :: k1 = 3./(1. + 2.*psi**(-0.5))    ,&
                       k2 = 10**(1.84148*(-log(psi))**0.5743)

    real :: d1,d2,d3,d4,d5


    DO j = ja,jz
       DO i = ia,iz
          DO k=1,m1  

             !- several particle/environment properties

             !- mean speed of air molecules (m/s)
             !  v_air = sqrt( 8. * kB	* temp3d(k,i,j) / (pi * M_AVEG) )
             v_air = SQRT( 7.3102e+2 * temp3d(k,i,j)  		      )

             !-dynamic viscosity of air (kg m^-1 s^-1)
             !  n_air = 1.8325e-5*(416.16/(temp3d(k,i,j)+120.))*(temp3d(k,i,j)/296.16)**1.5
             !optimized version
             n_air = 1.8325e-5*(416.16/(temp3d(k,i,j)+120.))*(temp3d(k,i,j)/296.16) * & 
                  SQRT(temp3d(k,i,j)/296.16)   
             !- kinematic viscosity of air ()
             nu = n_air/air_dens3d(k,i,j)
             !- mean free path of an air molecule (m)
             mfp = 2.* n_air /(air_dens3d(k,i,j)*v_air)

             DO ispc=naer_a,naer_z

                !print*,'particle=',ispc


                !- Knudsen number
                Kn = mfp/part_radius(ind_mode(ispc),ind_aer(ispc))

                !- Slip correction factor (Gi)
                Gi = 1. + Kn*( ASP + BSP*EXP(-CSP/Kn) )


                !- This is regime 1 (Pruppacher and Klett (chap. 10)) , first guess for terminal velocity
                !- 0.5 < part_radius < 10 micrometers  / 1.e-6 < Re < 1e-2
                !  
                !- particle sedimentation velocity (m/s)
                !  v_sed(k,i,j) = (2./9.)*g*part_radius**2  *(part_dens-air_dens3d(k,i,j))*Gi/n_air
                !- opt
                dd_sedim(ispc)%v_sed_part(k,i,j) = 2.18*part_radius(ind_mode(ispc),ind_aer(ispc))**2.  &
                     *(part_dens(ind_mode(ispc),ind_aer(ispc)))*Gi/n_air !part_dens >>dens_air


                ! Reynolds number
                !   re = 2*part_radius*v_sed_part(i,j)/nu
                re = 2. * air_dens3d(k,i,j) * part_radius(ind_mode(ispc),ind_aer(ispc)) * &
                     dd_sedim(ispc)%v_sed_part(k,i,j)/n_air

                go to 343
                !-----
		!- srf addded new formulation 
		!- Approach of FALL 3d (Costa et al., 2005)
		!
		!- avoid division by zero in case  re = zero.
		re = re + eps
		!- drag coefficient
		Cd = 24./(re * k1 ) * (1. + 0.1118 *( re *( k1*k2 )**0.6567 )) &
		     + 0.4305 * k2 / ( 1. + 3305. / (re * k1 * k2 ))
                 
                !-- particle sedimentation velocity (m/s)
		dd_sedim(ispc)%v_sed_part(k,i,j) = (4.* g *part_dens(ind_mode(ispc),ind_aer(ispc)) * &
		                                    2.  *part_radius(ind_mode(ispc),ind_aer(ispc))) /&
		                                   (3.* Cd * air_dens3d(k,i,j))
		
                dd_sedim(ispc)%v_sed_part(k,i,j) = ( dd_sedim(ispc)%v_sed_part(k,i,j) ) ** 0.5
		
		cycle ! not using the corrections below
		
		!------
 343            continue			 
                IF( re .GE. 0.01 .AND. re .LE. 300.) THEN
                   !  This is "regime 2" in Pruppacher and Klett (chap. 10).
                   x = LOG(24.*re/Gi)
                   y = -3.18657 + x*(0.992696    - x*(.00153193 &
                        + x*(0.000987059 + x*(.000578878&
                        - x*(8.55176E-05 - x* 3.27815E-06 )))))
                   !if( y .lt. -675. ) y = -675.
                   !if( y .ge.  741. ) y =  741.
                   y=MIN(MAX(y,-675.),741.)

                   re = EXP(y)*Gi

                   dd_sedim(ispc)%v_sed_part(k,i,j) = re * n_air / &
                        (2.*part_radius(ind_mode(ispc),ind_aer(ispc))* air_dens3d(k,i,j))
                ENDIF

                IF( re > 300. )THEN
                   !  This is "regime 3" in Pruppacher and Klett (chap. 10).
		   ! - units mus be CGS
		   d1= air_dens3d(k,i,j) *1.e-3 ! g/cm�
                   d2= part_dens(ind_mode(ispc),ind_aer(ispc)) *1.e-3 ! g/cm�
                   d3= g *100. ! cm/s2
                   d4=n_air*10. ! g/cm/s
                   ! - old
		   !z  = ((1.e6*air_dens3d(k,i,j)**2) &
                   !     /  (g * part_dens(ind_mode(ispc),ind_aer(ispc)) * n_air**4) )**(onesix)
 
                   z  = ((1.e6*d1**2)  /  (d3 * d2* d4**4) )**(onesix)

                   d5=dd_sedim(ispc)%v_sed_part(k,i,j)*100. !cm/s

                   !b0 = (24.*dd_sedim(ispc)%v_sed_part(k,i,j) *n_air)/100.
                   b0  = (24.*d5*d4)/100.
		   
		   x  = LOG(z*b0)
                   y  = -5.00015 + x*(5.23778   - x*(2.04914 - x*(0.475294 &
                        - x*(0.0542819 - x* 0.00238449 ))))
                   !if( y .lt. -675. )  y = -675.0
                   !if( y .ge.  741. )  y =  741.0
                   y=MIN(MAX(y,-675.),741.)
                   re = z*EXP(y)*Gi

                   dd_sedim(ispc)%v_sed_part(k,i,j) = re * n_air / &
                        (2.*part_radius(ind_mode(ispc),ind_aer(ispc)) * air_dens3d(k,i,j))

                ENDIF
                !tmp
                !for testing 
                !dd_sedim(ispc,ngrid)%v_sed_part(k,i,j) = float(ispc)/500.
                !if(k==1 .or. K==2)print*,'v-sed=',k,dd_sedim(ispc,ngrid)%v_sed_part(k,i,j) 
                !tmp	


             ENDDO ! ispc loop
          ENDDO  !   k loop    
       ENDDO
    ENDDO


  !print*,'max val min val '
  !do  ispc=naer_a,naer_z
  ! 
  ! if(ispc>7)   print*,1000.*part_radius(ind_mode(ispc),ind_aer(ispc)),&
  !	     maxval(dd_sedim(ispc)%v_sed_part),minval(dd_sedim(ispc)%v_sed_part)
  !    
  !enddo      
  !stop 333
    
      RETURN
  END SUBROUTINE sedim_particles_3d
  !------------------------------------------------------------------------
  SUBROUTINE fa_preptc_with_sedim(m1,m2,m3,vt3da,vt3db,vt3dc,vt3df    &
                                 ,vt3dk,dn0,rtgt,f13t,f23t,dtlt       &
                                 ,N_current_scalar,num_scalar_aer_1st &
                                 ,wp,wc,vt3dp,nzpmax,hw4,dzm,dzt      &
                                 ,dd_sedim)

    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3

    REAL    , INTENT(IN)    :: vt3da(m1,m2,m3)
    REAL    , INTENT(IN)    :: vt3db(m1,m2,m3)
    REAL    , INTENT(OUT)   :: vt3dc(m1,m2,m3)

    REAL    , INTENT(INOUT) :: vt3df(m1,m2,m3)

    REAL    , INTENT(INOUT) :: vt3dk(m1,m2,m3)
    REAL    , INTENT(IN)    :: dn0(m1,m2,m3)

    REAL    , INTENT(IN)    :: rtgt(m2,m3)

    REAL    , INTENT(IN)    :: f13t(m2,m3)
    REAL    , INTENT(IN)    :: f23t(m2,m3)

    !srf- aerosol section 
    REAL    , INTENT(IN)    :: dtlt
    INTEGER , INTENT(IN)    :: N_current_scalar
    INTEGER , INTENT(IN)    :: num_scalar_aer_1st

    REAL    , INTENT(IN)    :: wp(m1,m2,m3)
    REAL    , INTENT(IN)    :: wc(m1,m2,m3)
    REAL    , INTENT(INOUT) :: vt3dp(m1,m2,m3)

    ! grid_dims
    INTEGER , INTENT(IN)    :: nzpmax

    ! mem_grid
    REAL    , INTENT(IN)    :: hw4(nzpmax)
    REAL    , INTENT(IN)    :: dzm(nzpmax)
    REAL    , INTENT(IN)    :: dzt(nzpmax)

    TYPE(sedim_type), INTENT(INOUT) :: dd_sedim(naer_transported)

    INTEGER :: j,i,k,im,ip,jm,jp
    REAL    :: c1,c2,c3,c4,rtgti

    ! VT3DA, VT3DB, and VT3DC are input as the velocity components (averaged
    ! between past and current time levels) times dtlt.

    ! Add contribution to VT3DC from horiz winds crossing sloping sigma surfaces,
    !    and include 1/rtgt factor in VT3DC
    ! Compute half Courant numbers: VT3DD, VT3DE, and VT3DF
    ! Compute weight at scalar point: VT3DH
    ! Compute advective weights for the linear term: VT3DI, VCTR1, and VCTR2

    INTEGER ISPC

    !------------------ aerosol mapping
    ! num_scalar_aer_1st                  corresponds do aerosol  1 (=naer_a)
    ! num_scalar_aer_1st +1               corresponds do aerosol  2
    ! ...
    ! num_scalar_aer_1st +NAER_TRANSPORTED corresponds do aerosol naer_z
    ! in this case, to access the correct V_SED use "ISPC"
    ISPC = N_current_scalar - num_scalar_aer_1st + 1

    ! sedim is included at vertical direction (vt3dc)
    DO j = 1,m3
       DO i = 1,m2
          DO k = 1,m1
             !if(j==20 .and. i==20) print*,'particle',ISPC,sedim(ispc,ngrid)%v_sed_part(k,i,j),wp(k,i,j)
             vt3dc(k,i,j) = ( ( wp(k,i,j) + wc(k,i,j) )* 0.5 - dd_sedim(ispc)%v_sed_part(k,i,j) )* dtlt
             !-test	 vt3dc(k,i,j) = ( - sedim(ispc,ngrid)%v_sed_part(k,i,j) )* dtlt
          ENDDO
       ENDDO
    ENDDO
    !print*,'W=',wp(1:2,int(m2/2),int(m3/2)),wc(1:2,int(m2/2),int(m3/2))
!!!! (be sure the lines below are executed each timestep)!!!!
    !- only necessary one time 
    IF(ISPC==1) THEN
       DO j = 1,m3
          jm = MAX(1,j-1)
          !jp = min(m3,j+1)
          DO i = 1,m2
             im = MAX(1,i-1)
             !ip = min(m2,i+1)
             !rtgti = 1. / rtgt(i,j)
             DO k = 1,m1-1
                vt3dp(k,i,j) = ((vt3da(k,i,j) + vt3da(k+1,i,j)  &
                     + vt3da(k,im,j) + vt3da(k+1,im,j)) * f13t(i,j)  &
                     + (vt3db(k,i,j) + vt3db(k+1,i,j) + vt3db(k,i,jm)  &
                     + vt3db(k+1,i,jm)) * f23t(i,j)) * hw4(k)
                vt3dk(k,i,j) = dzt(k) / dn0(k,i,j)
             ENDDO

          ENDDO
       ENDDO
    ENDIF


    !-  then the fluxes are re-evaluated
    DO j = 1,m3
       jm = MAX(1,j-1)
       !jp = min(m3,j+1)
       DO i = 1,m2
          im = MAX(1,i-1)
          !ip = min(m2,i+1)
          rtgti = 1. / rtgt(i,j)
          DO k = 1,m1-1
             !-orig way
             !        vt3dc(k,i,j) = ((vt3da(k,i,j) + vt3da(k+1,i,j)  &
             !           + vt3da(k,im,j) + vt3da(k+1,im,j)) * f13t(i,j)  &
             !           + (vt3db(k,i,j) + vt3db(k+1,i,j) + vt3db(k,i,jm)  &
             !           + vt3db(k+1,i,jm)) * f23t(i,j)) * hw4(k)  &
             !           + vt3dc(k,i,j) * rtgti
             !opt way
             vt3dc(k,i,j) = vt3dp(k,i,j) &
                  + vt3dc(k,i,j) * rtgti


             vt3df(k,i,j) = .5 * vt3dc(k,i,j) * dzm(k)
             !vt3dk(k,i,j) = dzt(k) / dn0(k,i,j)

          ENDDO

       ENDDO
    ENDDO
    !do k = 1,m1-1
    !   vctr1(k) = (zt(k+1) - zm(k)) * dzm(k)
    !   vctr2(k) =  (zm(k) - zt(k)) * dzm(k)
    !enddo
    !print*,'W=',wp(1:2,int(m2/2),int(m3/2)),wc(1:2,int(m2/2),int(m3/2))

    !convert velocity components * dtlt (VT3DA, VT3DB, VT3DC)
    ! into mass fluxes times dtlt.

    DO j = 1,m3
       DO i = 1,m2
          !c1 = fmapui(i,j) * rtgu(i,j)
          !c2 = fmapvi(i,j) * rtgv(i,j)
          DO k = 1,m1-1
             vt3dc(k,i,j) = vt3dc(k,i,j) * .5  &
                  * (dn0(k,i,j) + dn0(k+1,i,j))      !air_dens3d ????
          ENDDO
       ENDDO
    ENDDO
  END SUBROUTINE fa_preptc_with_sedim
  !--------------------------------------------------------------------------



END MODULE module_dry_dep
