subroutine emiss(igrd,m1,m2,m3,i0,j0,ia,iz,ja,jz, & nspc,nspecem,nlu,ndpspc,densfac,lmap,lprobe,nsrc, & idsrc,isrc,jsrc,ncol,nrow,nlay,nlay_ems,deltat, & dx,dy,mapscl,height,depth,windu,windv,tempk,press, & icdocn,fsurf,tsurf,eflxnh3, & aremis,pttrace,numpts,armass, & ptmass,reemis,conc,depfld,ipa_cel) use bndary use chmstry use filunit use ptemiss use procan use tracer c c----CAMx v7Beta6 190902 c c EMISS updates species concentrations due to emissions c c Copyright 1996 - 2019 c Ramboll c c Modifications: c 10/29/01 Map scale factor added to emissions increment c 04/21/04 Added optional plume rise override c 04/08/10 Updates to plume rise algorithm (multi-layer injection) c 01/15/13 Allow user-set multi-layer injection using EFFPH and FLOWRAT c 04/29/13 Inject surface model re-emissions c 12/18/15 Add oceanic I2/HOI emissions from O3 dep c 07/23/18 Updated for Bi-Di NH3 drydep c c Input arguments: c igrd grid number c nspc number of species c nlu number of land use categories c ndpspc number of deposition species inc. I2,HOI emissions c densfac gas concentration conversion factor c lmap emissions source species map c lprobe flag for a Probing Tool call c nsrc number of point sources c idsrc point source ID map c isrc grid column index for point sources c jsrc grid row index for point sources c ncol number of columns c nrow number of rows c nlay number of layers c nlay_ems number of layers in emissions files c deltat time step size (s) c dx cell size in x-direction (m) c dy cell size in y-direction (m) c mapscl map scale factor c height layer height (m) c depth layer depth (m) c windu wind speed in x-direction (m/s) c windv wind speed in y-direction (m/s) c tempk air temprature (K) c press air pressure (mb) c icdocn ocean map c fsurf landuse fraction c tsurf surface temperature (K) c eflxnh3 emiss flux of Bi-Di NH3 drydep (umol/m2-s) c aremis area source strength (mol/s) c pttrace point source strength (mol/s) c numpts number of point sources c armass mass from area source emission (umol) c ptmass mass from point source emission (umol) c reemis emissions from surface model (mol/s) c conc species concentrations (umol/m3) c depfld deposited mass carrying I2,HOI emissions (mol) c ipa_cel gridded array to identify if cell is c in a IPRM sub-domain c c Output arguments: c conc species concentrations (umol/m3) c armass mass from area source emission (umol) c ptmass mass from point source emission (umol) c c Routines called: c PLUMERIS c IXEMIS c c Called by: c EMISTRNS c implicit none include "camx.prm" include "flags.inc" c c-----Input arguments c integer :: igrd integer :: m1,m2,m3,i0,j0,ia,iz,ja,jz integer :: nspc,nspecem,nsrc,ncol,nrow,nlay,nlay_ems,numpts integer :: nlu integer :: ndpspc real :: deltat,dy real :: densfac logical :: lprobe integer icdocn(m1,m2) integer lmap(*) real*8 armass(*) real*8 ptmass(*) integer idsrc(*) integer isrc(*) integer jsrc(*) real pttrace(numpts,nspecem) real conc(m1,m2,m3,nspc) real aremis(m1,m2,nlay_ems,nspecem) real reemis(m1,m2,nsmspc) real dx(nrow) real mapscl(m1,m2) real, dimension(m1,m2,m3) :: windu,windv,tempk,press real, dimension(m1,m2,m3) :: height, depth real, dimension(m1,m2) :: tsurf real, dimension(m1,m2) :: eflxnh3 real, dimension(m1,m2,nlu) :: fsurf real, dimension(m1,m2,ndpspc) :: depfld c c======================== Process Analysis Begin ==================================== c integer ipa_cel(ncol,nrow,nlay) integer ipa_idx c c========================= Process Analysis End ===================================== c c-----Local variables c integer l,i,j,k,kstk,kt,kb,lsrc,n integer ll,lsm,imod,iflgbidi,imap,ipart real frac real hght1d(0:MXLAYER) real wind1d(MXLAYER) real tempk1d(MXLAYER) real dtdz1d(MXLAYER) real pfrac(MXLAYER) real gamma,p0,grav real i2flx,hoiflx real cellarea,vol,dconc,w2,dz,dtheta,zstk,zbot,ztop,dstkabs,wp,tp,zrise, & trise,pwidth,rip,fp,qp2,cq1,cq2,rkp,bot,top,fwtr,co3,convfac real*8 dmass c data gamma,p0,grav /0.286,1000.,9.8/ data cq1,cq2 /0.4,3.0/ c c-----Entry point c c-----Update concentration due to area source c if (larsrc) then c c$omp parallel default(shared) c$omp& private(l,i,j,imod,ll,lsm,cellarea,vol,dmass,dconc,ipa_idx, c$omp& i2flx,hoiflx,fwtr,co3,convfac,w2,iflgbidi,imap,ipart,frac) c c$omp do schedule(static) c do l = 1,nspc imod = lmap(l) if( imod .LE. 0 .AND. l .NE. ki2 .AND. l .NE. khoi ) cycle if (lsrfmod) then lsm = 0 do ll = 1,nsmspc if (l .EQ. idsmsp(ll)) lsm = ll enddo endif iflgbidi = 0 if (lbidinh3) then if (lprobe) then !bk - assume that ITRNH3 is the only tracer group that includes NH3 (for SA) ! and NAMALL is the only species group that includes NH3 (for DDM/HDDM); ! this achieves better efficiency at the cost of some flexibility if (ltrace .AND. tectyp.NE.RTRAC .AND. tectyp.NE.RTCMC) then do ll = 1,ntrcls if (idxcls(ll).eq.ITRNH3) then do imap = 1,nregin if (l.eq.(iemcls(ll)-1+imap+(IGRPBIDI-1)*nregin)) & iflgbidi = imap enddo endif enddo else if (lddm .OR. lhddm) then do ll = 1,nemddm if (emddmsp(ll).eq.'NH3 ' .or. emddmsp(ll).eq.NAMALL) then do imap = 1,nregin if (l.eq.(iptddm(knh3)+nbdic+(ll-1)*ngroup*nregin+(imap-1)*ngroup+IGRPBIDI-1)) & iflgbidi = imap enddo endif enddo endif else if (l.eq.knh3) iflgbidi = 1 endif endif do i = 2, m1-1 do j = 2, m2-1 cellarea = dx(j+j0)*dy/(mapscl(i,j)**2.) vol = cellarea*depth(i,j,1) dmass = 0. do k=1,nlayers_ems c c------Standard gridded emissions c if( imod .GT. 0 ) dmass = DBLE(aremis(i,j,k,imod)*deltat*1.e6) c c------Re-emissions from Bi-Di NH3 drydep c if( iflgbidi .GT. 0 ) then if (lprobe) then if (ltrace .AND. tectyp.NE.RTRAC .AND. tectyp.NE.RTCMC) then do ipart = 1,npartial(IGRPBIDI,igrd) imap = igrmap(IGRPBIDI,ipart,igrd,i+i0,j+j0) frac = frcmap(IGRPBIDI,ipart,igrd,i+i0,j+j0) if (imap.EQ.iflgbidi) dmass = dmass & + DBLE(eflxnh3(i,j)*cellarea*deltat*frac) ! [umol] enddo else if (lddm .OR. lhddm) then imap = igrmap(0,1,igrd,i+i0,j+j0) if (imap.EQ.iflgbidi) dmass = dmass & + DBLE(eflxnh3(i,j)*cellarea*deltat) ! [umol] endif else dmass = dmass + DBLE(eflxnh3(i,j)*cellarea*deltat) ! [umol] endif endif c c------Re-emissions from surface model c if (lsrfmod .and. lsm .NE. 0) & dmass = dmass + DBLE(reemis(i,j,lsm)*deltat*1.d6) c c------Oceanic halogen emissions related to ozone deposition c if (k .eq. 1 .AND. lixemis .AND. icdocn(i,j).GT.0 & .AND. .NOT.lprobe .AND. & (l.EQ.ki2 .OR. l.EQ.khoi)) then w2 = sqrt(windu(i,j,1)**2. + windv(i,j,1)**2.) w2 = amax1(w2,0.1) if (nlu.EQ.NLUW89) then fwtr = fsurf(i,j,7) else fwtr = fsurf(i,j,1) endif convfac = densfac*(273./tempk(i,j,1))* & (press(i,j,1)/1013.) co3 = 1000.*conc(i,j,1,ko3)/convfac call ixemis(dx(j+j0),dy,mapscl(i,j),icdocn(i,j),fwtr, & co3,tsurf(i,j),w2,i2flx,hoiflx) if (l.eq.ki2) & dmass = dmass + DBLE(i2flx*deltat*1.d6) depfld(i,j,ndpspc-1) = depfld(i,j,ndpspc-1) + & i2flx*deltat/ & (dx(j+j0)*dy)*(mapscl(i,j)**2) if (l.eq.khoi) & dmass = dmass + DBLE(hoiflx*deltat*1.d6) depfld(i,j,ndpspc) = depfld(i,j,ndpspc) + & hoiflx*deltat/ & (dx(j+j0)*dy)*(mapscl(i,j)**2) endif dconc = REAL(dmass)/vol if( i .GE. ia .and. i .LE. iz .AND. & j .GE. ja .AND. j .LE. jz ) armass(l) = & armass(l) + dmass conc(i,j,k,l) = conc(i,j,k,l) + dconc enddo c c======================== Process Analysis Begin ==================================== c c----- if surface layer in this column is in a sub-domain c then track the area emissions in this grid cell --- c if( (.NOT. ltrace) .AND. lipr ) then if( i .GE. ia .and. i .LE. iz .AND. & j .GE. ja .AND. j .LE. jz ) then if( ipa_cel(i+i0,j+j0,1) .GT. 0 ) then ipa_idx = ipa_cel(i+i0,j+j0,1) cipr(IPR_AEMIS, ipa_idx, l) = & cipr(IPR_AEMIS, ipa_idx, l) + dconc endif endif endif c c========================= Process Analysis End ===================================== c enddo enddo enddo c c --- end of parallelized loop --- c c$omp end parallel c endif c c-----Update concentration due to point sources c if (lptsrc) then do 50 lsrc = 1,nsrc n = idsrc(lsrc) i = isrc(lsrc)-i0 j = jsrc(lsrc)-j0 if (i < 2 .or. i >= m1 .or. j < 2 .or. j >= m2) goto 50 c hght1d(0) = 0. do k = 1,nlay hght1d(k) = height(i,j,k) tempk1d(k) = tempk(i,j,k) w2 = windu(i,j,k)*windu(i,j,k) + windv(i,j,k)*windv(i,j,k) wind1d(k) = amax1(sqrt(w2),0.1) if (k.lt.nlay) then dz = height(i,j,k+1)/2. if (k.gt.1) dz = (height(i,j,k+1) - & height(i,j,k-1))/2. dtheta = (tempk(i,j,k+1)*(p0/press(i,j,k+1))**gamma - & tempk(i,j,k)*(p0/press(i,j,k))**gamma) dtdz1d(k) = dtheta/dz else dtdz1d(k) = dtdz1d(k-1) endif enddo c c-----Apply plume rise OVERRIDE if effph is negative c if (effph(n) .lt. 0. .and. flowrat(n) .ge. 0.) then zstk = abs(effph(n)) do kstk = 1,nlay if (hght1d(kstk).gt.zstk) goto 14 enddo kstk = nlay zstk = hght1d(kstk-1) + 1. 14 zbot = max(0.,zstk - 1.) ztop = min(hght1d(nlay),zstk + 1.) c c-----Handle unique situation for EPS3 point source files c else if(effph(n) .LT. 0. .AND. flowrat(n) .LT. 0. .AND. & effph(n) .EQ. flowrat(n) ) then zstk = abs(effph(n)) do kstk = 1,nlay if (hght1d(kstk).gt.zstk) goto 13 enddo kstk = nlay zstk = hght1d(kstk-1) + 1. 13 zbot = max(0.,zstk - 1.) ztop = min(hght1d(nlay),zstk + 1.) c c-----Apply plume layer distribution if stack height is also negative --- c else if(effph(n) .LT. 0. .AND. flowrat(n) .LT. 0. .AND. & effph(n) .NE. flowrat(n) ) then zstk = 0.5*( ABS( effph(n) ) + ABS( flowrat(n) ) ) do kstk = 1,nlay if (hght1d(kstk) .GT. zstk) goto 18 enddo kstk = nlay 18 zbot = MAX(0., MIN( ABS( flowrat(n) ), ABS( effph(n) ) ) ) ztop = MIN( hght1d(nlay), & MAX( ABS( effph(n) ), ABS( flowrat(n) ) ) ) else c c-----Calculate plume rise c dstkabs = abs(dstk(n)) call plumeris(nlay,hght1d,tempk1d,dtdz1d,wind1d,hstk(n), & dstkabs,tstk(n),vstk(n),zstk) do kstk = 1,nlay if (hght1d(kstk).gt.zstk) goto 15 enddo kstk = nlay zstk = hght1d(kstk-1) + 1. c c-----Determine plume depth c 15 wp = amax1(1.,vstk(n)/2.) tp = (tempk1d(kstk) + tstk(n))/2. zrise = zstk - hstk(n) trise = zrise/wp pwidth = sqrt(2.)*dstkabs rip = grav*pwidth*abs(tp - tempk1d(kstk))/ & (tempk1d(kstk)*wp*wp) fp = 1. + 4.*rip qp2 = fp*wp*wp*(cq1 + cq2*wind1d(kstk)*wind1d(kstk)/ & (wind1d(kstk)*wind1d(kstk) + wp*wp)) rkp = 0.15*pwidth*sqrt(qp2) pwidth = 3.*sqrt(pwidth*pwidth + 2.*rkp*trise) pwidth = max(1.,min(pwidth,zrise)) zbot = max(hght1d(0),zstk - pwidth/2.) ztop = min(hght1d(nlay),zstk + pwidth/2.) endif c c-----Calculate layers to receive emission injection c pwidth = ztop - zbot do k = 1,nlay pfrac(k) = 0. enddo do kt = kstk,nlay if (hght1d(kt).ge.ztop) goto 16 enddo kt = nlay 16 do kb = 1,kstk if (hght1d(kb).gt.zbot) goto 17 enddo 17 do k = kb,kt bot = max(hght1d(k-1),zbot) top = min(hght1d(k),ztop) pfrac(k) = (top - bot)/pwidth enddo do 40 l = 1,nspc imod = lmap(l) if( imod .LE. 0 ) goto 40 if (lpiglet(n) .AND. ipigflg .NE. 0) goto 40 do k = kb,kt vol = dx(j+j0)*dy*depth(i,j,k)/(mapscl(i,j)**2.) dmass = DBLE(pttrace(n,imod)*deltat*1e6*pfrac(k)) dconc = REAL(dmass)/vol if( i .GE. ia .and. i .LE. iz .AND. & j .GE. ja .AND. j .LE. jz ) ptmass(l) = & ptmass(l) + dmass conc(i,j,k,l) = conc(i,j,k,l) + dconc c c======================== Process Analysis Begin ==================================== c c --- if layer containing plume for this cell is in a sub-domain c then track the point source emissions for this cell ---- c if( (.NOT. ltrace) .AND. lipr ) then if( i .GE. ia .and. i .LE. iz .AND. & j .GE. ja .AND. j .LE. jz ) then if( ipa_cel(i+i0,j+j0,k) .GT. 0 ) then ipa_idx = ipa_cel(i+i0,j+j0,k) cipr(IPR_PTEMIS, ipa_idx, l) = & cipr(IPR_PTEMIS, ipa_idx, l) + dconc endif endif endif c c========================= Process Analysis End ===================================== c enddo 40 continue 50 continue endif c return end