      subroutine drydep(m1,m2,m3,itzon,
     &                  tsurf,cellat,cellon,lrdlai,lai,
     &                  pwc,cwc,height,press,windu,windv,cod,
     &                  water,fsurf,tempk,snow,icdocn,conc,vdep,eflxnh3)
      use filunit
      use chmstry
      use bndary
      use camxcom
c
      implicit none

c
c----CAMx v7Beta6 190902
c 
c     DRYDEP is the driver for the calculation of gridded dry deposition 
c     velocities for a given grid. Deposition velocities are calculated for
c     each gas species and for each aerosol size bin, weighted by the
c     fractional land use specified for each cell.
c 
c     Copyright 1996 - 2018
c     Ramboll
c           
c     Modifications:
c        4/4/00    Added aerosol deposition as f(size)
c        4/4/01    Fixed a few bugs in the call for VD_AER
c        1/9/02    Aerosol size and density now species-dependent
c        3/26/03   Added scaling factor to surface resistance (provided
c                  on chemparam file), and zero dry deposition for 
c                  insoluble gases
c        4/9/03    Removed rain, added precip and cloud water contents:
c                  surfaces can now be rain-wetted, fog-wetted, or dew-wetted
c        6/6/03    Protect against divide by zero with totland
c        6/11/03   Use optional surface roughness length, if available
c        7/21/03   Use optional drought stress and snow cover, if available;
c                  Introduced latitude-dependent specification of season;
c                  Revised solar flux calculation to use RADM cloud adjustment
c        8/18/03   Relate drought stress to Palmer Drought Index
c        4/21/04   Incorporated sectional PM
c        11/19/04  Incorporated season-dependent roughness length
c        05/11/09  Now checks that chemistry is on before accessing the
c                  kph2o pointer.
c        8/10/09   Added the Zhang et al (2003) dry deposition option
c        02/11/11  Removed optional roughness from AHO
c        04/02/12  Removed drought stress and snow flag; AHO
c                  file is now just ozone column
c        04/17/12  Special handling of Hg dep to ocean surface
c        08/11/14  Snow cover changed to water equivalent depth
c        12/07/14  Modified for VBS
c        07/23/18  Updated for Bi-Di NH3 drydep
c
c     Input arguments:
c        itzon               time zone
c        tsurf               surface temperature field (K)
c        cellat              cell centroid latitude (deg)
c        cellon              cell centroid longitude (deg)
c        lrdlai              flag that LAI data were read
c        lai                 gridded LAI data from file
c        pwc                 precipitation water content (g/m3)
c        cwc                 cloud water content (g/m3)
c        height              layer interface height field (m)
c        press               layer pressure field (mb)
c        windu               layer U-component wind field (m/s)
c        windv               layer V-component wind field (m/s)
c        cod                 cloud optical depth field (unitless)
c        water               layer water vapor field (ppm)
c        fsurf               fractional landuse cover field (fraction)
c        tempk               layer temperature field (K)
c        snow                snow cover water equivalent (m)
c        icdocn              gridded ocean mask
c        conc                concentration field (umol/m3, ug/m3)
c             
c     Output arguments: 
c        vdep                species-dependent deposition velocity field (m/s)
c        eflxnh3             emiss flux of Bi-Di NH3 drydep (umol/m2-s)
c             
c     Routines called: 
c        CALDATE
c        GETZNTH
c        MICROMET
c        VD_GAS
c        VD_AER
c        HENRYFNC
c             
c     Called by: 
c        EMISTRNS 
c 
      include 'camx.prm'
      include 'deposit.inc'
      include 'camx_aero.inc'
      include 'flags.inc'
c
      integer :: itzon, idate, month, i, j, k, m, mbin, latbin, isesn
      real    :: zenith, deltaz, temp0, prss0, ucomp, vcomp, wind
      real    :: qwatr, ev, es, rh, dew, totland, z0, ustar, el
      real    :: psih, wstar, henso2, henry, vd, qtf, qtc, vtf, vtc
      real    :: diamc, vdf, roprta, diam, vdc, kwtr, pbl, zl
      real    :: qt, vt
      real    :: convfac, conc_nh3, fnet
      integer :: iwet, l, iflgso2, iflgo3, iflgnh3, isemptyf, isemptyc
      integer :: isec, iaero, isempty
      integer :: m1,m2,m3,iday,iyear
      real, dimension(m1,m2) :: tsurf, cellat, cellon, lai, snow
      integer, dimension(m1,m2) :: icdocn
      integer, dimension(12) :: nday
      real, dimension(m1,m2,m3) :: height
      real, dimension(m1,m2,m3) :: press, windu, windv, cod,
     &                             water, tempk, pwc, cwc
      real, dimension(m1,m2,NLU) :: fsurf
      real, dimension(m1,m2,nspec) :: vdep
      real, dimension(m1,m2,m3,nspec) :: conc
      real, dimension(m1,m2) :: eflxnh3
      real, parameter :: eflxnh3_fac(NLUZ03) = (/      ! flag to shut off Bi-Di NH3 re-emissions
     &                                            1.0, !  water
     &                                            1.0, !  ice
     &                                            1.0, !  inland lake
     &                                            1.0, !  evergreen needleleaf trees
     &                                            1.0, !  evergreen broadleaf trees
     &                                            1.0, !  deciduous needleleaf trees
     &                                            1.0, !  deciduous broadleaf trees
     &                                            1.0, !  tropical broadleaf trees
     &                                            1.0, !  drought deciduous trees
     &                                            1.0, !  evergreen broadleaf shrub
     &                                            1.0, !  deciduous shrubs
     &                                            1.0, !  thorn shrubs
     &                                            1.0, !  short grass and forbs
     &                                            1.0, !  long grass
     &                                            0.0, !  crops
     &                                            0.0, !  rice
     &                                            0.0, !  sugar
     &                                            0.0, !  maize
     &                                            0.0, !  cotton
     &                                            0.0, !  irrigated crops
     &                                            1.0, !  urban
     &                                            1.0, !  tundra
     &                                            1.0, !  swamp
     &                                            1.0, !  desert
     &                                            1.0, !  mixed wood forests
     &                                            1.0  !  transitional forest
     &                                          /)

      real lai_ref_intpl,rlai,ref_lai,lai_f,pi
      real tau,ctrns,fcld
      real*8 coszen,solflux
      logical lstable
      logical ldark,lrdlai
c
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/
c
c-----Entry point
c
      pi = 3.1415927
      idate = date
      call caldate(idate)
      iyear = idate/10000
      month = (idate - 10000*iyear)/100
      if (mod(iyear,4).eq.0) nday(2)=29
      iday = idate - 100*(idate/100)
c
c-----Loop over rows and columns
c
c$omp parallel default(shared)
c$omp&  private(i,j,mbin,latbin,isesn,ref_lai,m,lai_ref_intpl,
c$omp&          rlai,zenith,ldark,coszen,tau,k,ctrns,fcld,
c$omp&          solflux,deltaz,temp0,prss0,
c$omp&          ucomp,vcomp,wind,qwatr,ev,es,rh,iwet,dew,
c$omp&          l,totland,z0,lai_f,ustar,el,psih,wstar,lstable,zl,
c$omp&          henso2,henry,iflgso2,iflgo3,iflgnh3,vd,isemptyf,isemptyc,qtf,
c$omp&          qtc,vtf,vtc,diamc,vdf,roprta,diam,vdc,kwtr,isec,qt,vt,
c$omp&          convfac,conc_nh3,fnet,
c$omp&          iaero,isempty)
c$omp do schedule(static)
c
      do 30 j = 2,m2-1 
        do 20 i = 2,m1-1
c
c-----Determine season
c
          mbin = month
          if (cellat(i,j).lt.0.) then
            mbin = mod(month+6,12)
            if (mbin.eq.0) mbin = 12
          endif
          latbin = 1
          if (abs(cellat(i,j)).gt.20.) then
            latbin = 2
          elseif (abs(cellat(i,j)).gt.35.) then
            latbin = 3
          elseif (abs(cellat(i,j)).gt.50.) then
            latbin = 4
          elseif (abs(cellat(i,j)).gt.75.) then
            latbin = 5
          endif
          if ((cellat(i,j).gt.50. .and. cellat(i,j).lt.75.) .and.
     &        (cellon(i,j).gt.-15. .and. cellon(i,j).lt.15.)) latbin = 3
          isesn = iseason(latbin,mbin)
c
c-----Use input snow cover to set season, if specified
c
          if (snow(i,j).ge.0.001) isesn = 4    ! Snow cover check > 1 cm
c
c-----Determine cell-average relative LAI
c
          if (lrdlai) then
            ref_lai = 0.
            do m = 1,NLU
              lai_ref_intpl = lai_ref(m,mbin) + 
     &                    float(min(nday(mbin),iday))/float(nday(mbin))*
     &                    (lai_ref(m,mbin+1) - lai_ref(m,mbin))
              ref_lai = ref_lai + fsurf(i,j,m)*lai_ref_intpl
            enddo
            rlai = lai(i,j)/(ref_lai + 1.e-10)
          endif
c
c-----Calculate solar flux
c
          call getznth(cellat(i,j),cellon(i,j),time,date,itzon,
     &                 zenith,ldark)
          coszen = cos(DBLE(zenith)*DBLE(pi)/180.)
          tau = cod(i,j,m3)
          do k = m3-1,1,-1
            tau = tau + cod(i,j,k)
          enddo
          if (tau.lt.5.) then
            ctrns = 1.
            fcld  = 0.
          else
            ctrns = (5. - exp(-tau))/(4. + 0.42*tau)
            fcld  = 1.
          endif
          ctrns = 1.6*ctrns*coszen

          solflux = (990.*coszen - 30.)*
     &              (1. - fcld*(1. - ctrns))
          solflux = dmax1(DBLE(0.),solflux)
c
c-----Load local met variables
c
          deltaz = height(i,j,1)/2.
          temp0 = tsurf(i,j) - 273.15
          prss0 = press(i,j,1) - 2.*deltaz*(press(i,j,2) - 
     &                            press(i,j,1))/height(i,j,2)
          ucomp = (windu(i,j,1) + windu(i-1,j,1))/2.
          vcomp = (windv(i,j,1) + windv(i,j-1,1))/2.
          wind = sqrt(ucomp**2. + vcomp**2.)
          wind = amax1(0.1,wind)
c
c-----Determine surface wetness
c
          qwatr = 1.e-6*water(i,j,1)*18./28.8
          ev = qwatr*prss0/(qwatr + eps) 
          es = e0*exp((lv/rv)*(1./273. - 1./tsurf(i,j))) 
          rh = amin1(1.,ev/es)
          iwet = 0
          if (pwc(i,j,1).gt.cwmin) then
            iwet = 2
          elseif (cwc(i,j,1).gt.cwmin) then
            iwet = 1
          else
            dew = (1. - rh)*(wind + 0.6) 
            if (dew.lt.0.19) iwet = 1
          endif
c
c-----Loop over land use; surface roughness for water is dependent on
c     wind speed
c
          do l = 1,nspec
            vdep(i,j,l) = 0.
          enddo
          totland = 0.
          if (lbidinh3) then
            convfac = densfac*(273./tempk(i,j,1))*(press(i,j,1)/1013.)
            conc_nh3 = MAX( conc(i,j,1,knh3), bdnl(knh3)*convfac )*17.0 ! [ug/m3]
            eflxnh3(i,j) = 0.0
          endif
c
          do 10 m = 1,NLU
            if (fsurf(i,j,m).lt.0.01) goto 10
            totland = totland + fsurf(i,j,m)
c
c-----Set surface roughness for Wesely (1989) scheme
c
            if (idrydep.eq.1) then
              z0 = z0lu(m,isesn)
              if (m.eq.7) z0 = amax1(z0,2.0e-6*wind**2.5)
c
c-----Set surface roughness and LAI for Zhang (2003) scheme
c
            else
              lai_f = lai_ref(m,mbin) +
     &                float(min(nday(mbin),iday))/float(nday(mbin))*
     &                (lai_ref(m,mbin+1) - lai_ref(m,mbin))
              if (lrdlai) then
                lai_f = lai_f*rlai
                lai_f = amin1(lai_ref(m,15),lai_f)
                lai_f = amax1(lai_ref(m,14),lai_f)
              endif
              if (m.eq.1 .or. m.eq.3) then
                z0 = 2.0e-6*wind**2.5
              else 
                if (z02(m).gt.z01(m)) then
                  z0 = z01(m) + (lai_f - lai_ref(m,14))/
     &                          (lai_ref(m,15) - lai_ref(m,14))*
     &                          (z02(m) - z01(m))
                else
                  z0 = z01(m)
                endif
              endif
            endif
c
c-----Get surface layer micrometeorological parameters for this cell and
c     landuse type
c 
            call micromet(tempk(i,j,1),tsurf(i,j),press(i,j,1),prss0,
     &                    deltaz,wind,z0,pbl,ustar,el,psih,wstar,lstable)
            zl = deltaz/el
            if (zl.lt.0.) then
              zl = min(-5.,zl)
            else
              zl = max(5.,zl)
            endif 
c
c-----Loop over GAS species, and calculate deposition velocity for this cell,
c     landuse, and current species.
c     Use input drought stress code, if specified
c
            call henryfnc(0,henso20,tfactso2,tsurf(i,j),7.,1,1,1,henso2)
            do 40 l = nrad+1,ngas
              if (henry0(l).gt.1.e-6) then
                call henryfnc(l,henry0(l),tfact(l),tsurf(i,j),7.,knh3,
     &                        khno3,kso2,henry)
                iflgso2 = 0
                iflgo3 = 0
                iflgnh3 = 0
                if (l.eq.kso2) then
                  iflgso2 = 1
                  henry = henso2
                endif
                if (l.eq.ko3) iflgo3 = 1
                if (lbidinh3 .and. l.eq.knh3) iflgnh3 = 1
                if (idrydep .eq. 1) then   ! Use Wesely (1989) algorithm
                  call vd_gas(m,0,iwet,iflgso2,iflgo3,z0,
     &                     deltaz,psih,ustar,diffrat(l),henry,henso2,
     &                     f0(l),rscale(l),temp0,dstress(0),solflux,
     &                     rj(m,isesn),rlu(m,isesn),rac(m,isesn),
     &                     rlcs(m,isesn),rlco(m,isesn),rgss(m,isesn),
     &                     rgso(m,isesn),icdocn(i,j),vd)
                elseif (idrydep .eq. 2) then ! Use Zhang (2003) algorithm
                  call vd_gas_zhang(deltaz,zl,z0,ustar,tempk(i,j,1),
     &                      tsurf(i,j),solflux,rh,fcld,
     &                      pwc(i,j,1),coszen,m,snow(i,j),iflgo3,iflgnh3,conc_nh3,
     &                      henry,henso2,f0(l),diffrat(l),rscale(l),lai_f,vd,fnet)
                  if (iflgnh3.eq.1 .and. fnet.gt.0.0) 
     &                      eflxnh3(i,j) = eflxnh3(i,j)
     &                                   + (fnet/17.0)*eflxnh3_fac(m)*fsurf(i,j,m) ! [umol/m2-s]
                endif
              else
                vd = 0.
              endif
              if (spname(l).eq.'HG0       ' .AND. icdocn(i,j).eq.1 .AND.
     &            ((idrydep.eq.1 .and. m.eq.7) .OR.
     &             (idrydep.eq.2 .and. m.eq.1)))
     &          vd = 0.
              vdep(i,j,l) = vdep(i,j,l) + vd*fsurf(i,j,m)
 40         continue
c
c-----Loop over AEROSOL size bins, and calculate deposition velocity for
c     this cell and landuse
c
            if (naero .gt. 0) then
c
c-----Recalculate particle size for wet diameter
c
c
c-----CF 2-bin scheme
c
              if (lchem .AND. aeropt.eq.'CF') then
                isemptyf = 1
                isemptyc = 1
                qtf = 0. ! fine dry total mass
                qtc = 0. ! coarse dry total mass
                vtf = 0. ! fine dry total volume
                vtc = 0. ! coarse dry total volume
                do l = ngas+1,nspec
                  if (dcut(l,2).lt.(dcut(kph2o,2)+1.e-5)) then     !fine
                    if (l.ne.kph2o) then
                      if (conc(i,j,1,l).gt.bdnl(l)) isemptyf = 0
                      qtf = qtf + conc(i,j,1,l)
                      vtf = vtf + conc(i,j,1,l)/roprt(l)
                    endif
                  else                                             !coarse
                    if (conc(i,j,1,l).gt.bdnl(l)) isemptyc = 0
                    qtc = qtc + conc(i,j,1,l)
                    vtc = vtc + conc(i,j,1,l)/roprt(l)
                    diamc = sqrt(dcut(l,1)*dcut(l,2))
                  endif
                enddo
                vdf = 0.
                if (isemptyf.eq.0) then
                  roprta = (qtf + conc(i,j,1,kph2o)) / 
     &                     (vtf + conc(i,j,1,kph2o)/roprt(kph2o))
                  diam = sqrt(dcut(kph2o,1)*dcut(kph2o,2))
                  diam = 1.e-6*diam*(1. + 
     &                   conc(i,j,1,kph2o)/roprt(kph2o)/vtf)**0.33333
                  if (idrydep.eq.1) then     ! Use Wesely (1989) algorithm
                    call vd_aer(z0,deltaz,psih,ustar,diam,roprta,
     &                        tsurf(i,j),vdf)
                  elseif (idrydep.eq.2) then ! Use Zhang (2001) algorithm
                     call vd_aer_zhang(deltaz,zl,z0,ustar,diam,roprta,
     &                           tsurf(i,j),tempk(i,j,1),m,lai_f,vdf)
                  endif
                endif
                vdc = 0.
                if (isemptyc.eq.0) then
                  roprta = qtc/vtc
                  diam = 1.e-6*diamc
                  if (idrydep.eq.1) then     ! Use Wesely (1989) algorithm
                    call vd_aer(z0,deltaz,psih,ustar,diam,roprta,
     &                        tsurf(i,j),vdc)
                  elseif (idrydep.eq.2) then ! Use Zhang (2001) algorithm
                     call vd_aer_zhang(deltaz,zl,z0,ustar,diam,roprta,
     &                           tsurf(i,j),tempk(i,j,1),m,lai_f,vdc)
                  endif
                endif
                do l = ngas+1,nspec
                  if (dcut(l,2).lt.(dcut(kph2o,2)+1.e-5)) then
                    vdep(i,j,l) = vdep(i,j,l) + vdf*fsurf(i,j,m)
                  else
                    vdep(i,j,l) = vdep(i,j,l) + vdc*fsurf(i,j,m)
                  endif
                enddo
c
c-----CMU multi-section scheme
c
              elseif (lchem .AND. aeropt.eq.'CMU') then
                kwtr = (kph2o_1 - ngas)/nbin + 1
                if (nbin.eq.1) kwtr = kph2o_1 - ngas
                do isec = 1,nbin
                  isempty = 1
                  qt = 0. ! dry total mass
                  vt = 0. ! dry total volume
                  do iaero = 1,naero
                    if (iaero.ne.kwtr) then
                      l = ngas + (iaero - 1)*nbin + isec
                      if (conc(i,j,1,l).gt.bdnl(l)) isempty = 0
                      qt = qt + conc(i,j,1,l)
                      vt = vt + conc(i,j,1,l)/roprt(l)
                    endif
                  enddo
                  vd = 0.
                  if (isempty.eq.0) then
                    roprta = (qt + conc(i,j,1,kph2o_1-1+isec)) /
     &                       (vt + conc(i,j,1,kph2o_1-1+isec)/
     &                        roprt(kph2o_1))
                    diam = sqrt(dcut(ngas+isec,1)*dcut(ngas+isec,2))
                    diam = 1.e-6*diam*(1. + conc(i,j,1,kph2o_1-1+isec)/
     &                     roprt(kph2o_1)/vt)**0.33333
                    if (idrydep.eq.1) then     ! Use Wesely (1989) algorithm
                      call vd_aer(z0,deltaz,psih,ustar,diam,roprta,
     &                          tsurf(i,j),vd)
                    elseif (idrydep.eq.2) then ! Use Zhang (2001) algorithm
                     call vd_aer_zhang(deltaz,zl,z0,ustar,diam,roprta,
     &                           tsurf(i,j),tempk(i,j,1),m,lai_f,vd)
                    endif
                  endif
                  do iaero = 1,naero
                    l = ngas + (iaero - 1)*nbin + isec
                    vdep(i,j,l) = vdep(i,j,l) + vd*fsurf(i,j,m)
                  enddo
                enddo
c
c-----Other scheme
c
              else
                do l = ngas+1,nspec
                  diam = 1.e-6*sqrt(dcut(l,1)*dcut(l,2))
                  if (idrydep.eq.1) then     ! Use Wesely (1989) algorithm
                    call vd_aer(z0,deltaz,psih,ustar,diam,roprt(l),
     &                        tsurf(i,j),vd)
                  elseif (idrydep.eq.2) then ! Use Zhang (2001) algorithm
                     call vd_aer_zhang(deltaz,zl,z0,ustar,diam,roprt(l),
     &                           tsurf(i,j),tempk(i,j,1),m,lai_f,vd)
                  endif
                  vdep(i,j,l) = vdep(i,j,l) + vd*fsurf(i,j,m)
                enddo
              endif
            endif
 10       continue
c
c-----Calculate final landuse-weighted deposition velocities
c
          totland = amax1(totland, 0.01)
          do l = 1,nspec
            vdep(i,j,l) = vdep(i,j,l)/totland
          enddo
          if (lbidinh3) eflxnh3(i,j) = eflxnh3(i,j)/totland
c
 20     continue
 30   continue
c
c$omp end parallel
c
      return
      end
