subroutine drydeprt(m1,m2,m3,nrtsp,itzon, & tsurf,cellat,cellon,lrdlai,lai,pwc,cwc,height, & press,windu,windv,cod,water,fsurf,tempk, & snow,icdocn,vdep) use filunit use chmstry use bndary use camxcom use rtracchm use tracer c c----CAMx v7Beta6 190902 c c This version is for the RTRAC species. c DRYDEPRT 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 7/22/03 Taken from DRYDEP for regular model c 11/19/04 Incorporated season-dependent roughness length c 7/19/10 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 08/11/14 Snow cover changed to water equivalent depth c 07/23/18 Updated VD_GAS_ZHANG call arguments for Bi-Di NH3 c c Input arguments: c ncol number of columns c nrow number of rows c nlay number of layers c nrtsp number of species 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 c Output arguments: c vdep species-dependent deposition velocity field (m/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 :: m1,m2,m3,iday,iyear integer, dimension(12) :: nday real, dimension(m1,m2) :: tsurf, cellat, cellon, lai, snow integer, dimension(m1,m2) :: icdocn 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 lai_ref_intpl,rlai,ref_lai,lai_f real tau,ctrns,fcld real*8 coszen,solflux logical ldark,lrdlai,lstable integer :: iflgnh3 real :: rdum1, rdum2 c data pi/3.1415927/ data nday/31,28,31,30,31,30,31,31,30,31,30,31/ c c-----Entry point c 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 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 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.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,nrtrac vdep(i,j,l) = 0. enddo totland = 0. 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 call henryfnc(0,henso20,tfactso2,tsurf(i,j),7.,1,1,1, & henso2) do 40 l = 1,nrtgas if (rthlaw(l).gt.1.e-6) then call henryfnc(0,rthlaw(l),rttfact(l),tsurf(i,j),7.,1, & 1,1,henry) iflgso2 = 0 iflgo3 = 0 iflgnh3 = 0 if (idrydep .eq. 1) then ! Use Wesely (1989) algorithm call vd_gas(m,0,iwet,iflgso2,iflgo3,z0, & deltaz,psih,ustar,rtdrate(l),henry,henso2, & rtreact(l),rtscale(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,rdum1,henry,henso2,rtreact(l), & rtdrate(l),rtscale(l),lai_f,vd,rdum2) endif else vd = 0. endif 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 (nrtaero .gt. 0) then do 50 l = nrtgas+1,nrtrac diam = sqrt(rtlcut(l)*rtucut(l))*1.e-6 if (idrydep.eq.1) then ! Use Wesely (1989) algorithm call vd_aer(z0,deltaz,psih,ustar,diam,rtdens(l), & tsurf(i,j),vd) elseif (idrydep.eq.2) then ! Use Zhang (2001) algorithm call vd_aer_zhang(deltaz,zl,z0,ustar,diam,rtdens(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) 50 continue endif 10 continue c c-----Calculate final landuse-weighted deposition velocities c totland = amax1(totland, 0.01) do l = 1,nrtrac vdep(i,j,l) = vdep(i,j,l)/totland enddo c 20 continue 30 continue c return end