subroutine wetdep(m1,m2,m3,i0,j0,ixcl,jycl,igrd,ncol,nrow,nlay, & nspcs,deltat,area,depth,tempk,press, & cwc,pwr,pws,pwg,cph,densfac,dtout,ipa_cel,conc, & fluxtmp,depfld1,depfld2,wetfld,iptrsa) use chmstry use procan use tracer c c----CAMx v7Beta6 190902 c c WETDEP modifies vertical concentration profiles for a given grid via c precipitation processes. This subroutine has been completely rewritten c for CAMx v4. c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 06/15/03 Fixed bug in scaling of tmass c 02/11/04 Updated empirical rainfall relationships c 04/21/04 Incorporated sectional PM c 05/03/05 Revised gas scavenging calculation c 05/04/05 Revised DDM adjustments c 06/20/05 Revised to handle liquid/frozen cloud and precip water c 10/21/05 Added dissociation to Henry's Law for NH3, HNO3, SO2 c 08/31/06 Added map scale factor c 07/16/07 Fixed PA bug - adjust cipr if evap c 07/16/07 -bkoo- Added check for HDDM c 05/12/08 -cemery- Revised deposited liquid concentration to use c correct rainfall volume c 07/16/08 -bkoo- Added DDM turn-off flag c 07/16/08 -gwilson- Now accumulates mass flux summary outside OMP loop c 12/22/09 -cemery- Improved treatment of "c0" array c 09/23/13 -bkoo- Placed levap-setting code outside of the spc loop c Added DDM-PM code c 10/09/14 -cemery- Restructured for subgrid convective model c 12/07/14 -bkoo- Modified for VBS c 03/17/16 -cemery- Fixed bug to dissallow tmass<0 (Ondrej Vlcek, CHMI) c 06/06/16 -cemery- Simplified gas scavenging approach, removed c rain evap c 08/04/16 -cemery- Revised mean PM size in the fine mode (<2.5um) to c account for more efficient PM scavenging at the c high end of the size distribution c c Input arguments: c m1 number of slice columns c m2 number of slice rows c m3 number of slice layers c i0 slice i-offset c j0 slice j-offset c ixcl i-cell index c jycl j-cell index c igrd grid index c ncol number of grid columns c nrow number of grid rows c nlay number of grid layers c nspcs number of total species c deltat time step (s) c area cell area (m2) c depth cell depth (m) c tempk temperature field (K) c press pressure field (mb) c cwc cloud water content (g/m3) c pwr rain water content (g/m3) c pws snow water content (g/m3) c pwg graupel water content (g/m3) c cph cloud water pH c densfac factor to convert to umol/m3 c dtout output frequency (minutes) c ipa_cel identify if cell is in a IPRM sub-domain c conc concentration field (umol/m3, ug/m3) c c Output arguments: c conc concentration field (umol/m3, ug/m3) c fluxtmp temporary array for fluxes c depfld1 wet deposited mass (mol/ha, g/ha) c depfld2 surface liquid concentrations (mol/l, g/l) c wetfld wet depostion array for tracer species c c Routines called: c SCAVRAT c HENRYFNC c c Called by: c CIGDRIVE c implicit none include 'camx.prm' include 'flags.inc' c integer nlay,nspcs,igrd real deltat,area,densfac,dtout real, dimension(nlay) :: tempk,press,cwc,pwr,pws,pwg,cph,depth real, dimension(nspcs) :: depfld1,depfld2,fluxtmp real, dimension(MXLAYER,MXSPEC) :: conc c logical lcloud,ltop,lgraupl integer kbot,ktop,k,kzcl,l,isemptyf,isemptyc,kwtr, & isec,isempty,iaero real rd,rhoh2o,cellvol,rhoair,delc,delm,convfac,cmin,hlaw, & dscav,gscav,ascav,qtf,qtc,vtf,vtc,psizec,ascavf, & roprta,psize,ascavc,qt,vt,cwat,pwat,pwtr,rconst,volume, & dtfall,drpvel,delc0,ceqc,ceqp,ceqpp,pliq c real pp real rr real volrat real c0(MXSPEC) real tmass(MXSPEC) real delr(MXSPEC) c c======================== Probing Tool Begin =========================== c integer m1,m2,m3,i0,j0,ixcl,jycl integer ncol,nrow integer icls integer*8 iptrsa real fc2r,fr2c,fc2rc0,fr2cc0,cyctr logical lzerc c real c0trac(MXTRSP) real tmtrac(MXTRSP) real delcls(MXALCLS) real c0cls(MXALCLS) real concls(MXALCLS) real cytcls(MXALCLS) c integer ipa_cel(ncol,nrow,nlay), ipa_idx real, dimension(m1,m2,notimespc) :: wetfld c c========================= Probing Tool End ============================ c data rd /287./ ! Dry air gas constant (J/K/kg) data rhoh2o /1.e6/ ! water density (g/m3) data rconst /8.206e-2/ ! gas constant (l.atm/mol.K) c c-----Entry point c c-----If precip is not getting to the ground, there is no wet dep c if (pwr(1).lt.cwmin .and. pws(1).lt.cwmin .and. & pwg(1).lt.cwmin) return c c-----Scan column for layers containing cloud top c kbot = 1 pp = pwr(1) + pws(1) + pwg(1) volrat = pp/rhoh2o ! drop volume/air volume rr = (volrat/1.0e-7)**1.27 ! rainfall rate (mm/hr) c ktop = 0 do k = 1,nlay if (cwc(k).lt.cwmin .and. ktop.ne.0) then goto 26 endif if (cwc(k).ge.cwmin) ktop = k enddo ktop = nlay 26 continue c do l = 1,nspec c0(l) = 0. tmass(l) = 0. enddo c c-----Loop over layers and species for precipitating columns c ltop = .true. pliq = 0. do 30 k = ktop,kbot,-1 kzcl = k lcloud = .false. lgraupl = .false. if (cwc(k).ge.cwmin) lcloud = .true. if (pwg(k).ge.cwmin) lgraupl = .true. cellvol = area*depth(k) rhoair = 100.*press(k)/(rd*tempk(k)) cwat = cwc(k) pwtr = pp + pliq if (tempk(k).lt.273.) & cwat = max(0.,cwat*(tempk(k) - tamin)/(273. - tamin)) c c======================== Probing Tool Begin =========================== c if( ((lddm .OR. lhddm) .AND. lddmcalc(igrd)) .OR. & (ltrace .AND. tectyp.NE.RTRAC .AND. tectyp.NE.RTCMC) ) then if (ltop) then do l = 1,ntotsp tmtrac(l) = 0. c0trac(l) = 0. enddo endif do icls = 1,ntrcls delcls(icls) = 0. concls(icls) = 0. c0cls(icls) = 0. cytcls(icls) = 0. enddo endif c c======================== Probing Tool End ============================= c c-----Calculate scavenging for soluble gas species c do 40 l = nrad+1,ngas if( henry0(l).LT.1.e-6 ) goto 40 convfac = densfac*(273./tempk(k))*(press(k)/1013.) cmin = bdnl(l)*convfac conc(k,l) = max(cmin,conc(k,l)) call henryfnc(l,henry0(l),tfact(l),tempk(k),cph(k),knh3, & khno3,kso2,hlaw) hlaw = hlaw*rconst*tempk(k) ceqc = conc(k,l)*(1. - 1./(1. + hlaw*cwat/rhoh2o)) call henryfnc(l,henry0(l),tfact(l),tempk(k),5.,knh3, & khno3,kso2,hlaw) hlaw = hlaw*rconst*tempk(k) pwat = pwtr if (tempk(k).lt.273. .and. rscale(l).gt.0.) pwat = 0. ceqp = conc(k,l)*(1. - 1./(1. + hlaw*pwat/rhoh2o)) ceqpp = (conc(k,l) - ceqc)*(1. - 1./(1. + hlaw*pwat/rhoh2o)) c call scavrat(.false.,lcloud,lgraupl,tamin,rr,tempk(k), & rhoair,diffrat(l),0.,0.,dscav,gscav,ascav,drpvel) c dtfall = depth(k)/drpvel if (c0(l).lt.ceqp .or. ceqp.eq.0.) then delc = ceqc*(1. - exp(-dscav*deltat)) + & ceqpp*(1. - exp(-gscav*deltat)) delc0 = ceqc*(1. - exp(-dscav*min(dtfall,deltat))) + & ceqpp*(1. - exp(-gscav*min(dtfall,deltat))) if (ceqp.eq.0.) then delc = min(delc,conc(k,l)-cmin) delc0 = min(delc0,conc(k,l)-cmin) else delc = min(delc,ceqp-c0(l),conc(k,l)-cmin) delc0 = min(delc0,ceqp-c0(l),conc(k,l)-cmin) endif else delc = (ceqp - c0(l))*(1. - exp(-gscav*deltat)) delc0 = (ceqp - c0(l))*(1. - exp(-gscav*min(dtfall,deltat))) endif delm = delc*cellvol if (delm.lt.0.) then delm = max(delm,-tmass(l)) delc = delm/cellvol endif c c======================== Probing Tool Begin =========================== c c --- Accumulate OSAT/PSAT tracer adjustments by tracer class --- c if( ltrace .AND. tectyp.NE.RTRAC .AND. tectyp.NE.RTCMC ) then cyctr = 0.5*min(conc(k,l),c0(l)) if (pwat.eq.0) then cyctr = cyctr*(1. - exp(-dscav*deltat)) else cyctr = cyctr*(1. - exp(-(dscav+gscav)*deltat)) endif do icls = 1,ntrcls concls(icls) = concls(icls) + conc(k,l)*fluxmap(l,icls) c0cls(icls) = c0cls(icls) + c0(l)*fluxmap(l,icls) delcls(icls) = delcls(icls) + delc*fluxmap(l,icls) cytcls(icls) = cytcls(icls) + cyctr*fluxmap(l,icls) enddo endif c c --- For DDM convert delc to flux --- c if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then if( delc.GT.0. ) then fc2r = delc/conc(k,l) fr2c = 0. fc2rc0 = delc0/conc(k,l) fr2cc0 = 0. elseif( c0(l) .GT. 1.0e-20 ) then fc2r = 0. fr2c = -delc/c0(l) fc2rc0 = 0. fr2cc0 = -delc0/c0(l) else fc2r = 0. fr2c = 0. fc2rc0 = 0. fr2cc0 = 0. endif lzerc = .FALSE. if( conc(k,l)-delc .LE. 1.1*cmin ) lzerc = .TRUE. call adjddmc0(l,ixcl,jycl,kzcl,fc2r,fr2c,fc2rc0,fr2cc0, & c0trac,tmtrac,cellvol,lzerc, & m1,m2,m3,ntotsp,ptconc(iptrsa)) endif c c-----PA change from wet deposition c if( lipr ) then if( ipa_cel(ixcl+i0,jycl+j0,k) .GT. 0 ) then ipa_idx = ipa_cel(ixcl+i0,jycl+j0,k) cipr(IPR_WDEP,ipa_idx,l) = cipr(IPR_WDEP,ipa_idx,l) - delc endif endif c c========================= Probing Tool End ============================ c c-----Update the cell concentration and rain mass c conc(k,l) = conc(k,l) - delc c0(l) = max(c0(l)+delc0,0.) tmass(l) = max(tmass(l)+delm,0.) 40 continue pliq = pliq + cwat*(1. - exp(-dscav*min(dtfall,deltat))) c c-----Calculate scavenging for particulate species c if (naero .gt. 0) then c c-----Recalculate particle size for wet diameter c Adjust geometric mean particle diameters <2.5um by 10x to account for c mass scavenging over the whole size distribution (assuming geometric c standard deviation of 2.0) 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(k,l).gt.bdnl(l)) isemptyf = 0 qtf = qtf + conc(k,l) vtf = vtf + conc(k,l)/roprt(l) endif else ! coarse if (conc(k,l).gt.bdnl(l)) isemptyc = 0 qtc = qtc + conc(k,l) vtc = vtc + conc(k,l)/roprt(l) psizec = sqrt(dcut(l,1)*dcut(l,2)) endif enddo ascavf = 0. if (isemptyf.eq.0) then roprta = (qtf + conc(k,kph2o)) / & (vtf + conc(k,kph2o)/roprt(kph2o)) psize = sqrt(dcut(kph2o,1)*dcut(kph2o,2)) psize = 1.e-6*psize*(1. + & conc(k,kph2o)/roprt(kph2o)/vtf)**0.33333 if (psize.lt.2.5e-6) psize = psize*10. call scavrat(.true.,lcloud,lgraupl,tamin,rr,tempk(k), & rhoair,0.,psize,roprta,dscav,gscav,ascavf,drpvel) endif ascavc = 0. if (isemptyc.eq.0) then roprta = qtc/vtc psize = psizec*1.e-6 call scavrat(.true.,lcloud,lgraupl,tamin,rr,tempk(k), & rhoair,0.,psize,roprta,dscav,gscav,ascavc,drpvel) endif do l = ngas+1,nspec if (dcut(l,2).lt.(dcut(kph2o,2)+1.e-5)) then delr(l) = 1. - exp(-ascavf*deltat) else delr(l) = 1. - exp(-ascavc*deltat) 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(k,l).gt.bdnl(l)) isempty = 0 qt = qt + conc(k,l) vt = vt + conc(k,l)/roprt(l) endif enddo ascav = 0. if (isempty.eq.0) then roprta = (qt + conc(k,kph2o_1-1+isec))/ & (vt + conc(k,kph2o_1-1+isec)/ & roprt(kph2o_1)) psize = sqrt(dcut(ngas+isec,1)*dcut(ngas+isec,2)) psize = 1.e-6*psize*(1. + conc(k,kph2o_1-1+isec)/ & roprt(kph2o_1)/vt)**0.33333 call scavrat(.true.,lcloud,lgraupl,tamin,rr,tempk(k), & rhoair,0.,psize,roprta,dscav,gscav,ascav,drpvel) endif do iaero = 1,naero l = ngas + (iaero-1)*nbin + isec delr(l) = 1. - exp(-ascav*deltat) enddo enddo c c-----Other scheme c else do l = ngas+1,nspec psize = 1.e-6*sqrt(dcut(l,1)*dcut(l,2)) if (psize.lt.2.5e-6) psize = psize*10. call scavrat(.true.,lcloud,lgraupl,tamin,rr,tempk(k), & rhoair,0.,psize,roprt(l),dscav,gscav,ascav,drpvel) delr(l) = 1. - exp(-ascav*deltat) enddo endif c do 50 l = ngas+1,nspec delc = 0. delm = 0. cmin = bdnl(l) conc(k,l) = max(cmin,conc(k,l)) delc = conc(k,l)*delr(l) delc = min(delc,conc(k,l)-cmin) c conc(k,l) = conc(k,l) - delc delm = delc*cellvol tmass(l) = tmass(l) + delm c c========================= Probing Tool Begin ============================ c c-----PA change from wet deposition c if( lipr ) then if( ipa_cel(ixcl+i0,jycl+j0,k) .GT. 0 ) then ipa_idx = ipa_cel(ixcl+i0,jycl+j0,k) cipr(IPR_WDEP,ipa_idx,l) = cipr(IPR_WDEP,ipa_idx,l) - & delc endif endif c if( ltrace .AND. tectyp.NE.RTRAC .AND. & tectyp.NE.RTCMC ) then do icls = 1,ntrcls delcls(icls) = delcls(icls) + delc*fluxmap(l,icls) enddo endif c c --- Adjust DDM sensitivities for PM species --- c if ( lddm .AND. lddmcalc(igrd) ) then fc2r = delr(l) fr2c = 0. fc2rc0 = 0. fr2cc0 = 0. lzerc = .FALSE. if( conc(k,l) .LE. 1.1*cmin ) lzerc = .TRUE. call adjddmc0(l,ixcl,jycl,kzcl,fc2r,fr2c,fc2rc0,fr2cc0, & c0trac,tmtrac,cellvol,lzerc, & m1,m2,m3,ntotsp,ptconc(iptrsa)) endif c c========================= Probing Tool End ============================ c 50 continue endif c c======================== Source Apportion Begin ======================= c c --- Adjust OSAT/PSAT tracer by tracer class --- c if( ltrace .AND. tectyp.NE.RTRAC .AND. tectyp.NE.RTCMC ) then call adjstc0(m1,m2,m3,igrd,ixcl,jycl,kzcl,delcls,concls, & cytcls,c0cls,tmtrac,cellvol) endif c c======================== Source Apportion End ========================= c ltop = .false. 30 continue c c-----Increment deposition flux arrays c volume = 3.6e-6*rr*deltat*area do l = nrad+1,nspec fluxtmp(l) = -tmass(l) depfld1(l) = 1.e-2*tmass(l)/area depfld2(l) = 1.e-9*(tmass(l)/volume)*deltat/(60.*dtout) enddo c c======================== Source Apportion Begin ======================= c if( lptdepout) then do l = 1,notimespc wetfld(ixcl,jycl,l) = wetfld(ixcl,jycl,l) + & 1.e-2*tmtrac(l)/area enddo endif c c======================== Source Apportion End ======================= c return end