subroutine cigdrive(m1,m2,m3,i0,j0,ia,iz,ja,jz,igrd,ncol,nrow, & nlay,nspcs,nspdep,idfin,deltat,densfac, & dtout,deltax,deltay,mapscl,ldark,depth,tempk, & press,water,cigfrc,cigwtr,cigph,cigpcp,cwc, & cph,pwr,pws,pwg,conc,cigmas,depfld,fluxes, & xmschem,iptrsa,ipa_cel,wetfld,fsurf,rtconc, & solmass,vegmass,m1rt,m2rt,m3rt,nrtsp) use chmstry use filunit use tracer use rtracchm c c----CAMx v7Beta6 190902 c c CIGDRIVE is the transport and chemistry driver for the cloud-in-grid c algorithm. It splits vertical column concentration profiles into in-cloud c and ambient, applies convective mass transport matrices to both, calls c aqueous chemistry and wet deposition for both, then re-combines final c profiles. c c If CiG is not invoked, this routine just performs wet deposition for c the regular gridded concentrations. c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 10/13/17 -bkoo- Added a dummy argument to AEROCHEM_CF call (aerosol pH) c 10/16/17 -bkoo- Added a dummy argument to AEROCHEM_CF call (time-weighted JNO2) c Changed DUMARR2 size 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 ia slice computational starting i-index c iz slice computational ending i-index c ja slice computational starting j-index c jz slice computational ending j-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 nspdep number of deposition species c idfin map of nested grids in this grid c deltat time step (s) c densfac factor to convert umol/m3 to ppm c dtout output frequency (minutes) c deltax cell size in x-direction (m) c deltay cell size in y-direction (m) c mapscl map scale factor c ldark darkness flag (T=dark) c depth cell depth (m) c tempk temperature field (K) c press pressure field (mb) c water humidity field (ppm) c cigfrc CiG column fractional coverage c cigwtr CiG cloud water (g/m3) c cigph CiG cloud pH c cigpcp CiG precipitation water (g/m3) c cph cloud water pH 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 conc concentration field (umol/m3, ug/m3) c cigmas CiG tracer transport matrix c c Output arguments: c conc concentration field (umol/m3, ug/m3) c depfld 2-D array of wet deposited mass (mol/ha, g/ha) c and surface liquid concentrations (mol/l, g/l) c fluxes boundary mass fluxes (umol, ug) c xmschem change in grid mass from chemistry (umol) c c Routines called: c SCAVRAT c HENRYFNC c c Called by: c EMISTRNS c implicit none include 'camx.prm' include 'flags.inc' include 'deposit.inc' include 'soap.inc' include 'chmdbg.inc' c integer m1,m2,m3,i0,j0,ia,iz,ja,jz integer igrd,ncol,nrow,nlay,nspcs,nspdep integer idfin(m1,m2) real deltat,densfac,dtout real deltax(nrow),deltay,mapscl(m1,m2) real depth(m1,m2,m3) real tempk(m1,m2,m3) real press(m1,m2,m3) real water(m1,m2,m3) real cigfrc(m1,m2) real cigwtr(m1,m2,m3) real cigph(m1,m2,m3) real cigpcp(m1,m2,m3) real cwc(m1,m2,m3) real cph(m1,m2,m3) real pwr(m1,m2,m3) real pws(m1,m2,m3) real pwg(m1,m2,m3) real conc(m1,m2,m3,nspcs) real depfld(m1,m2,3*nspdep) real cigmas(2,m3,m3,m1,m2) real*8 fluxes(nspcs,11),xmschem(nspcs) logical ldark(m1,m2) integer i,j,k,l,kk,ll,ks,ke,n,n1,n2 integer aero_flag real dtchem,area,acig,aamb,tcell,pcell,wcell, & convfac,cldph,cliq,volume,vol real tmassi,tmassf,terr real ftmp(MXSPEC),depfld1(MXSPEC),depfld2(MXSPEC) real con(MXLAYER,MXSPEC,2) real cnc(MXLAYER,MXSPEC) real dfac(MXLAYER) real cc(MXSPEC) real dtmp(MXSPEC) real dz(MXLAYER),tt(MXLAYER),pp(MXLAYER),cwtr(MXLAYER), & ph(MXLAYER),rwtr(MXLAYER),swtr(MXLAYER),gwtr(MXLAYER) real dumvar real dumarr1(NNV),dumarr2(7,MXSPEC),dumarr3(MXFDDM*NNV),dumarr4(1) real fluxtmp(MXCELLS,MXCELLS,MXSPEC) real dchem(MXCELLS,MXCELLS,MXSPEC) common /comcigdrive/ fluxtmp, dchem c c======================== Probing Tool Begin =========================== c integer*8 iptrsa integer nrtsp,m1rt,m2rt,m3rt integer ipa_cel(ncol,nrow,nlay) real wetfld(m1,m2,notimespc) real rtconc(m1rt,m2rt,m3rt,nrtsp) real fsurf(m1,m2,nlu) real solmass(m1rt,m2rt,nrtsp),vegmass(m1rt,m2rt,nrtsp) real cncrt(MXLAYER,MXTRSP) c c========================= Probing Tool End ============================ c c-----Entry point c do l = 1,nspec do j = 1,m2 do i = 1,m1 fluxtmp(i,j,l) = 0. dchem(i,j,l) = 0. enddo enddo enddo dtchem = deltat/3600. igrdchm = igrd aero_flag = 0 if (aeropt.EQ.'CF') aero_flag = 1 c c-----Loop over rows and columns c do 10 j = ja,jz !2,nrow-1 c c$omp parallel default(shared) c$omp& private(i,l,ll,n,n1,n2,ftmp,depfld1,depfld2,k,kk,ks,ke,con,area, c$omp& dz,tt,pp,acig,cnc,cwtr,ph,rwtr,swtr,gwtr,aamb,cncrt,tcell,pcell, c$omp& wcell,convfac,cc,cldph,cliq,tmassi,tmassf,terr,dfac,volume,vol, c$omp& dtmp,dumvar,dumarr1,dumarr2,dumarr3,dumarr4) c$omp& copyin(/ijkgrd/) c c$omp do schedule(dynamic) c do 20 i = ia,iz !2,ncol-1 if (idfin(i,j).gt.igrd) goto 20 do n = 1,2 do l = 1,nspec do k = 1,nlay con(k,l,n) = 0. enddo enddo enddo do k = 1,nlay dfac(k) = press(i,j,k)/tempk(i,j,k) enddo c c-----Apply convective transfer matrices to split the column profile c if (lcig(igrd) .AND. cigfrc(i,j).gt.0.) then do l = nrad+1,nspec tmassi = 0. tmassf = 0. do ks = 1,nlay !Starting layer tmassi = tmassi + conc(i,j,ks,l)*depth(i,j,ks) do ke = 1,nlay !Ending layer con(ke,l,1) = con(ke,l,1) + & conc(i,j,ks,l)*cigmas(1,ke,ks,i,j)*dfac(ke)/dfac(ks) con(ke,l,2) = con(ke,l,2) + & conc(i,j,ks,l)*cigmas(2,ke,ks,i,j)*dfac(ke)/dfac(ks) enddo enddo c c-----Check for conc < 0 and mass conservation c do k = 1,nlay if (con(k,l,1).lt.0. .or. con(k,l,2).lt.0.) then write(iout,'(//,a)') 'ERROR in CIGDRIVE:' write(iout,'(2a)') 'Negative concentration after', & ' convective transport' write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k write(iout,*) 'Spec Name In-cloud Ambient' do kk = 1,nlay write(iout,'(i3,2x,a7,2e10.3)') kk,spname(l), & con(kk,l,1),con(kk,l,2) enddo call camxerr() endif tmassf = tmassf + depth(i,j,k)*(cigfrc(i,j)*con(k,l,1) & + (1. - cigfrc(i,j))*con(k,l,2)) enddo terr = abs(tmassf - tmassi)/tmassi if (terr.gt.0.002) then write(iout,'(//,a)')'ERROR in CIGDRIVE:' write(iout,*)'CiG transport matrix not conserving mass' write(iout,*)'Grid,I,J : ',igrd,i,j write(iout,*)'Species : ',spname(l) write(iout,*)'Starting/ending mass: ',tmassi,tmassf call camxerr() endif enddo else do l = nrad+1,nspec do k = 1,nlay con(k,l,2) = conc(i,j,k,l) enddo enddo endif c c-----Wet deposition c if (lwet) then area = deltax(j+j0)*deltay/mapscl(i,j)**2 do k = 1,nlay dz(k) = depth(i,j,k) tt(k) = tempk(i,j,k) pp(k) = press(i,j,k) enddo c c-----CIG profile first c if (lcig(igrd) .AND. cigfrc(i,j).gt.0.) then acig = cigfrc(i,j)*area do l = nrad+1,nspec ftmp(l) = 0. depfld1(l) = 0. depfld2(l) = 0. do k = 1,nlay cnc(k,l) = con(k,l,1) enddo enddo do k = 1,nlay cwtr(k) = cigwtr(i,j,k) ph(k) = cigph(i,j,k) rwtr(k) = 0. swtr(k) = 0. gwtr(k) = 0. if (tt(k).ge.273.) then rwtr(k) = cigpcp(i,j,k) else gwtr(k) = cigpcp(i,j,k) endif enddo call wetdep(m1,m2,m3,i0,j0,i,j,igrd,ncol,nrow,nlay,nspec, & deltat,acig,dz,tt,pp,cwtr,rwtr,swtr,gwtr,ph, & densfac,dtout,ipa_cel,cnc,ftmp, & depfld1,depfld2,wetfld,iptrsa) c c-----Accumulate core model deposition fluxes c do l = nrad+1,nspec do k = 1,nlay con(k,l,1) = cnc(k,l) enddo fluxtmp(i,j,l) = fluxtmp(i,j,l) + ftmp(l) do ll = 1,ndepspc if (l .eq. ldepmap(ll)) then depfld(i,j,ndepspc+ll) = depfld(i,j,ndepspc+ll) + & depfld1(l) depfld(i,j,2*ndepspc+ll) = depfld(i,j,2*ndepspc+ll) & + depfld2(l) goto 100 endif enddo 100 continue enddo c endif c c-----Ambient profile next c aamb = area if (lcig(igrd) .AND. cigfrc(i,j).gt.0.) & aamb = (1. - cigfrc(i,j))*area do l = nrad+1,nspec ftmp(l) = 0. depfld1(l) = 0. depfld2(l) = 0. do k = 1,nlay cnc(k,l) = con(k,l,2) enddo enddo do k = 1,nlay cwtr(k) = cwc(i,j,k) ph(k) = cph(i,j,k) rwtr(k) = pwr(i,j,k) swtr(k) = pws(i,j,k) gwtr(k) = pwg(i,j,k) enddo call wetdep(m1,m2,m3,i0,j0,i,j,igrd,ncol,nrow,nlay,nspec, & deltat,aamb,dz,tt,pp,cwtr,rwtr,swtr,gwtr,ph, & densfac,dtout,ipa_cel,cnc,ftmp, & depfld1,depfld2,wetfld,iptrsa) c c-----Accumulate core model deposition fluxes c do l = nrad+1,nspec do k = 1,nlay con(k,l,2) = cnc(k,l) enddo fluxtmp(i,j,l) = fluxtmp(i,j,l) + ftmp(l) do ll = 1,ndepspc if (l .eq. ldepmap(ll)) then depfld(i,j,ndepspc+ll) = depfld(i,j,ndepspc+ll) + & depfld1(l) depfld(i,j,2*ndepspc+ll) = depfld(i,j,2*ndepspc+ll) + & depfld2(l) goto 101 endif enddo 101 continue enddo c c======================== Probing Tool Begin =========================== c c-----Wetdep for RTRAC c if( ltrace .AND. & (tectyp.EQ.RTRAC .OR. tectyp.EQ.RTCMC)) then do l = 1,nrtrac do k = 1,nlay cncrt(k,l) = rtconc(i,j,k,l) enddo enddo c call wetdeprt(m1,m2,m1rt,m2rt,i,j,nlay,nrtrac,deltat,aamb, & dz,tt,pp,cwtr,rwtr,swtr,gwtr,ph,densfac, & cncrt,solmass,vegmass,fsurf) c do l = 1,nrtrac do k = 1,nlay rtconc(i,j,k,l) = cncrt(k,l) enddo enddo endif endif c c========================= Probing Tool End ============================ c c-----Aqueous chemistry c if (lcig(igrd) .and. lchem .and. aeropt.EQ.'CF') then n1 = 2 n2 = 2 if (cigfrc(i,j).gt.0.) n1 = 1 do k = 1,nlay ichm = i jchm = j kchm = k c tcell = tempk(i,j,k) pcell = press(i,j,k) wcell = water(i,j,k) tchm = tcell wchm = wcell convfac = densfac*(273./tcell)*(pcell/1013.) volume = deltax(j+j0)*deltay*depth(i,j,k)/(mapscl(i,j)**2) c do n = n1,n2 !Loop over in-cloud and ambient profiles if (n.eq.1) then cldph = cigph(i,j,k) cliq = cigwtr(i,j,k) vol = cigfrc(i,j)*volume else cldph = cph(i,j,k) cliq = cwc(i,j,k) vol = volume if (cigfrc(i,j).gt.0.) vol = (1. - cigfrc(i,j))*volume endif if (tcell.lt.273.) & cliq = amax1(0.,cliq*(tcell - tamin)/(273. - tamin)) c c-----Sum mass and convert conc units before chemsitry c do l = 1,nspec dtmp(l) = con(k,l,n)*vol enddo do l = nrad+1,ngas cc(l) = con(k,l,n)/convfac enddo do l = ngas+1,nspec cc(l) = con(k,l,n) enddo do l = nrad+1,nspec if (cc(l).lt.0.) then write(iout,'(//,a)') 'ERROR in CIGDRIVE:' write(iout,*) 'Negative concentration before chem' write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k do ll = nrad+1,nspec write(iout,'(i3,2x,a7,e10.3)')ll,spname(ll),cc(ll) enddo call camxerr() endif cc(l) = amax1(bdnl(l),cc(l)) enddo c call aerochem_cf(wcell,tcell,pcell,cliq,cldph, & dumvar,dumvar, & dumarr1,cc,convfac,dtchem,dtchem, & dumarr2,dumarr3,dumarr4,1,1,.FALSE.,1, & aero_flag,.FALSE.,1,ldark(i,j)) c c-----Update cloud pH, convert conc units, and sum mass after chemsitry c if (n.eq.1) then cigph(i,j,k) = cldph else cph(i,j,k) = cldph endif do l = nrad+1,ngas con(k,l,n) = amax1(cc(l)*convfac,bdnl(l)) enddo do l = ngas+1,nspec con(k,l,n) = amax1(cc(l),bdnl(l)) enddo do l = 1,nspec dtmp(l) = con(k,l,n)*vol - dtmp(l) dchem(i,j,l) = dchem(i,j,l) + dtmp(l) enddo c enddo enddo endif c c-----Update main concentration array c if (lcig(igrd) .AND. cigfrc(i,j).gt.0.) then do l = nrad+1,nspec do k = 1,nlay conc(i,j,k,l) = cigfrc(i,j)*con(k,l,1) + & (1. - cigfrc(i,j))*con(k,l,2) enddo enddo else do l = nrad+1,nspec do k = 1,nlay conc(i,j,k,l) = con(k,l,2) enddo enddo endif c 20 continue c c --- end of parallelized loop --- c c$omp end parallel c 10 continue c c --- load fluxes into global array --- c do l = 1,nspec do j = ja,jz do i = ia,iz !2,ncol-1 fluxes(l,11) = fluxes(l,11) + fluxtmp(i,j,l) xmschem(l) = xmschem(l) + dchem(i,j,l) enddo enddo enddo c call flush(6) c return end