c*** OSATSA c subroutine osatsa(numcol,numrow,numlay,igrid,icell,jcell,kcell, & prdO3N,prdO3V,desto3,nspc,delcon, & convfac,cold,rrxn_irr,dtime) use grid use chmstry use tracer implicit none c c----CAMx v7Beta6 190902 c c----------------------------------------------------------------------- c Description: c This routine makes the "chemistry" adjustments to the tracer c species. The adjustments are based on the differences in c concentrations of the regular model species before and after c the regular model chemistry. This is essentially an adjustment c for production or decay. c c Copyright 1996 - 2018 c Ramboll c c Argument descriptions: c Inputs: c numcol I number of columns in this slice c numrow I number of rows in this slice c numlay I number of layers in this slice c igrid I grid number c icell I the X grid location of current cell c jcell I the X grid location of current cell c kcell I the vertical grid location of current layer c prdO3N R ozone production attributed to NOx c prdO3V R ozone production attributed to VOC c desto3 R ozone destruction c nspc I number of model species c delcon R array of change in each species contrntrations c convfac R conversion factor used for lower bound value c cold R array of concentrations at last time step c rrxn_irr R array of reactions from last chemistry step c dtime R change in time for current time step c c----------------------------------------------------------------------- c LOG: c----------------------------------------------------------------------- c c 05/22/97 --gwilson-- Now re-calculates the sum of VOC after c adjustment using wtkoh. c 09/20/03 --gwilson-- Changed the individual species args to c a vector c 07/20/04 --gwilson-- Changed for OSAT2 (ozone destruction) c 07/27/15 --gyarwood- Changed for OSAT3 (NO2 apportioned) c 09/11/15 --bkoo-- Merged APCA c 08/25/16 --gyarwood- Updated for CB6r4 c 11/09/16 --cemery-- Added Baker APCA point source override option c 10/16/17 --bkoo-- Changed DELCON size c c----------------------------------------------------------------------- c Include files: c----------------------------------------------------------------------- c include 'camx.prm' c c----------------------------------------------------------------------- c Argument declarations: c----------------------------------------------------------------------- c integer numcol integer numrow integer numlay integer igrid integer icell integer jcell integer kcell real prdO3N real prdO3V real desto3 integer nspc real delcon(7,nspc+1) real convfac real cold(*) real rrxn_irr(*) real dtime c c----------------------------------------------------------------------- c Local variables: c----------------------------------------------------------------------- c integer*8 idxcel, idx integer*8 idx_NIT, idx_HN3, idx_TPN, idx_NTR, idx_OO integer i, i_N, iHN3, iTPN, iNTR, iOO, iVOC integer istr1,istr2 real ptc1d(MXTRSP) real sumNIT, sumHN3, sumTPN, sumNTR, sumRGN real sumOO, sumO3, sumkOH, sumMIR, delVOC real destRGN, prodRGN, NIT2RGN, RGN2NIT real alphaRGN, betaRGN, gammaRGN, alphaNIT, cycNIT real alphaNTR, betaNTR, alphaHN3 real alphaTPN, betaTPN, gammaTPN, cycTPN real pairNOX, prodOO, destOO, OOtoO3, pairOO real bioNIT, bioVOC, facbio, O3used c c----------------------------------------------------------------------- c Entry point: c----------------------------------------------------------------------- c c --- calculate the index of the cell in the grid --- c idxcel = ipsa3d(igrid)-1 + DBLE(icell) + DBLE(numcol) * & DBLE(jcell-1) + DBLE(numcol) * DBLE(numrow) * & DBLE(kcell-1) istr1 = 4 istr2 = 6 if (lapcapt) then istr1 = 7 istr2 = 9 endif do i = 1, ntotsp idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptc1d(i) = ptconc(idx) enddo c c --- get the sum of the tracer species for N ---- c sumNIT = 0. bioNIT = 0. do i=iptcls(idxipt(ITRNIT)),nptcls(idxipt(ITRNIT)) sumNIT = sumNIT + ptc1d(i) if( ptname(i)(istr1:istr2) .EQ. '001' ) bioNIT = bioNIT + ptc1d(i) enddo sumNIT = MAX(BNDLPT, sumNIT) bioNIT = MAX(BNDLPT, bioNIT) c c --- get the sum of the tracer species for HN3 ---- c sumHN3 = 0. do i=iptcls(idxipt(ITRHN3)),nptcls(idxipt(ITRHN3)) sumHN3 = sumHN3 + ptc1d(i) enddo sumHN3 = MAX(BNDLPT, sumHN3) c c --- get the sum of the tracer species for TPN ---- c sumTPN = 0. do i=iptcls(idxipt(ITRTPN)),nptcls(idxipt(ITRTPN)) sumTPN = sumTPN + ptc1d(i) enddo sumTPN = MAX(BNDLPT, sumTPN) c c --- get the sum of the tracer species for NTR ---- c sumNTR = 0. do i=iptcls(idxipt(ITRNTR)),nptcls(idxipt(ITRNTR)) sumNTR = sumNTR + ptc1d(i) enddo sumNTR = MAX(BNDLPT, sumNTR) c c --- get the sum of the tracer species for RGN ---- c sumRGN = 0. do i=iptcls(idxipt(ITRRGN)),nptcls(idxipt(ITRRGN)) sumRGN = sumRGN + ptc1d(i) enddo sumRGN = MAX(BNDLPT, sumRGN) c c --- call routine to calculate the cycle coefficients --- c call cyctpnsa(nspec,convfac,cold,delcon,rrxn_irr, & destRGN,prodRGN,NIT2RGN,RGN2NIT, & alphaRGN,betaRGN,gammaRGN,alphaNIT,cycNIT, & alphaNTR,betaNTR,alphaHN3,alphaTPN,betaTPN, & gammaTPN,cycTPN) c c --- adjust RGN, NIT, HN3, TPN and NTR tracers ---- c i_N = iptcls(idxipt(ITRNIT)) - 1 iHN3 = iptcls(idxipt(ITRHN3)) - 1 iTPN = iptcls(idxipt(ITRTPN)) - 1 iNTR = iptcls(idxipt(ITRNTR)) - 1 c do i=iptcls(idxipt(ITRRGN)),nptcls(idxipt(ITRRGN)) i_N = i_N + 1 iHN3 = iHN3 + 1 iTPN = iTPN + 1 iNTR = iNTR + 1 c pairNOX = ptc1d(i) + ptc1d(i_N) c idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptconc(idx) = & sumRGN * ( (1.0 - cycNIT) * ptc1d(i) / sumRGN + & (cycNIT * pairNOX / (sumNIT + sumRGN)) ) & + NIT2RGN * ptc1d(i_N) / sumNIT & - (RGN2NIT + alphaRGN + betaRGN + & gammaRGN) * ptc1d(i) / sumRGN & + alphaNTR * ptc1d(iNTR) / sumNTR & + alphaHN3 * ptc1d(iHN3) / sumHN3 & + betaTPN * ptc1d(iTPN) / sumTPN c idx_NIT = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i_N-1) ptconc(idx_NIT) = & sumNIT * ( (1.0-cycNIT) * ptc1d(i_N) / sumNIT + & (cycNIT * pairNOX / (sumNIT + sumRGN)) ) & + RGN2NIT * ptc1d(i) / sumRGN & - (NIT2RGN + alphaNIT) * ptc1d(i_N) / sumNIT c idx_HN3 = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(iHN3-1) ptconc(idx_HN3) = ptconc(idx_HN3) & + gammaRGN * ptc1d(i) / sumRGN & + betaNTR * ptc1d(iNTR) / sumNTR & + alphaTPN * ptc1d(iTPN) / sumTPN & - alphaHN3 * ptc1d(iHN3) / sumHN3 c idx_TPN = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(iTPN-1) ptconc(idx_TPN) = ptconc(idx_TPN) & + betaRGN * ptc1d(i) / sumRGN & - (alphaTPN + betaTPN + gammaTPN) * ptc1d(iTPN) / sumTPN c idx_NTR = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(iNTR-1) ptconc(idx_NTR) = ptconc(idx_NTR) & + alphaNIT * ptc1d(i_N) / sumNIT & + alphaRGN * ptc1d(i) / sumRGN & + gammaTPN * ptc1d(iTPN) / sumTPN & - (alphaNTR + betaNTR) * ptc1d(iNTR) / sumNTR c ptconc(idx) = MAX(BNDLPT,ptconc(idx)) ptconc(idx_NIT) = MAX(BNDLPT,ptconc(idx_NIT)) ptconc(idx_HN3) = MAX(BNDLPT,ptconc(idx_HN3)) ptconc(idx_TPN) = MAX(BNDLPT,ptconc(idx_TPN)) ptconc(idx_NTR) = MAX(BNDLPT,ptconc(idx_NTR)) enddo c c --- if no O3 tracer, return here --- c if( .NOT. lozone ) return c c --- get the sum of the tracer species for O3N/O3V ---- c sumO3 = 0. do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3V)) sumO3 = sumO3 + ptc1d(i) enddo sumO3 = MAX(BNDLPT, sumO3) c c --- get the sum of the tracer species for OON/OOV ---- c sumOO = 0. do i=iptcls(idxipt(ITROON)),nptcls(idxipt(ITROOV)) sumOO = sumOO + ptc1d(i) enddo sumOO = MAX(BNDLPT, sumOO) c c --- set the local variable for VOC species concs --- c delVOC = 0. do i=1,nspec if( lvocsp(i) ) delVOC = delVOC + delcon(1,i) * crbnum(i) enddo c c --- loop over the VOC tracer species, calculate reactivity --- c sumkOH = 0. sumMIR = 0. bioVOC = 0. do i=iptcls(idxipt(ITRVOC)),nptcls(idxipt(ITRVOC)) sumkOH = sumkOH + ptc1d(i) * wtkoh(i) sumMIR = sumMIR + ptc1d(i) * wtmir(i) if( ptname(i)(istr1:istr2) .EQ. '001' ) & bioVOC = bioVOC + ptc1d(i) * wtmir(i) enddo sumkOH = MAX(BNDLPT, sumkOH) sumMIR = MAX(BNDLPT, sumMIR) bioVOC = MAX(BNDLPT, bioVOC) c c ---- decay VOC tracers --- c do i=iptcls(idxipt(ITRVOC)),nptcls(idxipt(ITRVOC)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptconc(idx) = ptconc(idx) + & delVOC * ptc1d(i) * wtkoh(i) / sumkOH ptconc(idx) = MAX(BNDLPT,ptconc(idx)) enddo c c --- update OON/OOV and allocate any ozone from NO2 destruction --- c prodOO = prodRGN destOO = destRGN OOtoO3 = MIN( destOO, prdO3N + prdO3V ) if( prdO3V .GT. 0.) then prdO3V = prdO3V - OOtoO3 else prdO3N = prdO3N - OOtoO3 endif c iOO = iptcls(idxipt(ITROON)) - 1 do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3V)) iOO = iOO + 1 idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) idx_OO = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(iOO-1) pairOO = ptc1d(i) + ptc1d(iOO) ptconc(idx) = & sumO3 * ( (1.0 - cycNIT) * ptc1d(i) / sumO3 + & (cycNIT * pairOO / (sumO3 + sumOO)) ) & + OOtoO3 * ptc1d(iOO) / sumOO ptconc(idx_OO) = & sumOO * ( (1.0-cycNIT) * ptc1d(iOO) / sumOO + & (cycNIT * pairOO / (sumO3 + sumOO)) ) & - destOO * ptc1d(iOO) / sumOO & + prodOO * ptc1d(i) / sumO3 ptconc(idx_OO) = MAX(BNDLPT,ptconc(idx_OO)) enddo c c --- allocate any ozone destruction across all tracers --- c if( desto3 .LT. 0. ) then do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3V)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptconc(idx) = ptconc(idx) + desto3 * ptc1d(i) / sumO3 enddo endif c c --- APCA c if ( lapca ) then c c --- set the biogenics contribution factor --- c facbio = MIN( bioNIT/sumNIT, bioVOC/sumMIR ) c c --- allocate VOC sensitive ozone production based on MIRs, c increase the O3V tracers except where limited by APCA --- c if( prdO3V .GT. 0. ) then O3used = 0. iVOC = iptcls(idxipt(ITRVOC)) - 1 do i=iptcls(idxipt(ITRO3V)),nptcls(idxipt(ITRO3V)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) iVOC = iVOC + 1 c c --- add O3 to biogenics group (group 1) based on it's c mimimum potential contribution --- c if( ptname(i)(istr1:istr2) .EQ. '001' ) then ptconc(idx) = ptconc(idx) + facbio * prdO3V * & ptc1d(iVOC) * wtmir(iVOC) / bioVOC O3used = O3used + facbio * prdO3V * & ptc1d(iVOC) * wtmir(iVOC) / bioVOC c c --- add O3 to other groups based on contribution of VOC conc --- c else ptconc(idx) = ptconc(idx) + prdO3V * & ptc1d(iVOC) * wtmir(iVOC) / sumMIR O3used = O3used + prdO3V * & ptc1d(iVOC) * wtmir(iVOC) / sumMIR endif enddo c c --- add any ozone not accounted for to O3N tracers, based on c contribution to anthropognic NOx --- c if( sumNIT .GT. bioNIT .AND. prdO3V .GT. O3used ) then i_N = iptcls(idxipt(ITRNIT)) - 1 do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3N)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) i_N = i_N + 1 if( ptname(i)(istr1:istr2) .NE. '001' ) then ptconc(idx) = ptconc(idx) + (prdO3V - O3used) * & ptc1d(i_N) / (sumNIT - bioNIT) endif enddo endif endif c c --- allocate NOx sensitive ozone production, c increase the O3N tracers except where limited by APCA --- c if( prdO3N .GT. 0. ) then O3used = 0. i_N = iptcls(idxipt(ITRNIT)) - 1 do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3N)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) i_N = i_N + 1 c c --- add O3 to biogenics group (group 1) based on it's c mimimum potential contribution --- c if( ptname(i)(istr1:istr2) .EQ. '001' ) then ptconc(idx) = ptconc(idx) + facbio * prdO3N * & ptc1d(i_N) / bioNIT O3used = O3used + facbio * prdO3N * & ptc1d(i_N) / bioNIT c c --- add O3 to other groups based on contribution of NOx conc --- c else ptconc(idx) = ptconc(idx) + prdO3N * & ptc1d(i_N) / sumNIT O3used = O3used + prdO3N * ptc1d(i_N) / sumNIT endif enddo c c --- add any ozone not accounted for to O3V tracers, based on c contribution to anthropognic VOC --- c if( sumMIR .GT. bioVOC .AND. prdO3N .GT. O3used ) then iVOC = iptcls(idxipt(ITRVOC)) - 1 do i=iptcls(idxipt(ITRO3V)),nptcls(idxipt(ITRO3V)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) iVOC = iVOC + 1 if( ptname(i)(istr1:istr2) .NE. '001' ) then ptconc(idx) = ptconc(idx) + (prdO3N - O3used) * & ptc1d(iVOC) * wtmir(iVOC) / (sumMIR - bioVOC) endif enddo endif endif c c --- OSAT c else c c --- allocate VOC sensitive ozone production based on MIRs, c increase the O3V tracers --- c if( prdO3V .GT. 0. ) then iVOC = iptcls(idxipt(ITRVOC)) - 1 do i=iptcls(idxipt(ITRO3V)),nptcls(idxipt(ITRO3V)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) iVOC = iVOC + 1 ptconc(idx) = ptconc(idx) + prdO3V * & ptc1d(iVOC) * wtmir(iVOC) / sumMIR enddo endif c c --- allocate NOx sensitive ozone production, c increase the O3N tracers --- c if( prdO3N .GT. 0. ) then i_N = iptcls(idxipt(ITRNIT)) - 1 do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3N)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) i_N = i_N + 1 ptconc(idx) = ptconc(idx) + prdO3N * & ptc1d(i_N) / sumNIT enddo endif c endif ! APCA? c c --- lower bound of ozone tracers c do i=iptcls(idxipt(ITRO3N)),nptcls(idxipt(ITRO3V)) idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptconc(idx) = MAX(BNDLPT,ptconc(idx)) enddo c c --- decay the timing tracer species --- c if(ntrtim .GT. 0 ) then do i=ipttim+1,nsaspc,2 idx = idxcel + DBLE(numcol) * DBLE(numrow) * & DBLE(numlay) * DBLE(i-1) ptconc(idx) = ptconc(idx) * EXP( -0.08333333 * dtime ) ptconc(idx) = MAX(BNDLPT,ptconc(idx)) enddo endif c c----------------------------------------------------------------------- c Return point: c----------------------------------------------------------------------- c return end