subroutine soachem(iout,tempk,press,dthr, & avgox,cprec,cCG,cNV,yfac, & nfam,avgsen,sen_prec,sen_CG,sen_NV,lddm) c c----CAMx v7Beta6 190902 c c----------------------------------------------------------------------- c Description: c This routine solves the oxidation reactions for SOA precursors c assuming first-order decay and constant oxidant concentrations c c c Rate constant: c k = A exp(-E/T) ; A [cm3/molecule-sec] @ E [K] c c Copyright 1996 - 2018 c Ramboll c c Input arguments: c iout standard output file unit c tempk cell temperature [K] c press cell pressure [mb] c dthr time duration to be integrated [hr] c avgox average oxidant concentrations [ppm] c 1 - O; 2 - OH; 3 - O3; 4 - NO3; 5 - NO; 6 - HO2 c cprec precursor concentrations [ppm] c 1 - BNZA; 2 - TOLA; 3 - XYLA; 4 - IVOA c 5 - ISP ; 6 - TRP ; 7 - SQT c cCG CG species concentrations [ppm] c 1 - CG1; 2 - CG2; 3 - CG3; 4 - CG4 c cNV NV products concentration [ppm] c 1 - CGA (anthro); 2 - CGB (bio) c yfac weighted yield factors c nfam number of sensitivity families c avgsen average sensitivities of oxidants [ppm] c sen_prec sensitivities of precursors [ppm] c sen_CG sensitivities of CGs [ppm] c sen_NV sensitivities of NVs [ppm] c lddm logical flag for DDM sensitivities c c Output arguments: c cprec precursor concentrations [ppm] c cCG CG species concentrations [ppm] c cNV NV products concentration [ppm] c yfac weighted yield factors c sen_prec sensitivities of precursors [ppm] c sen_CG sensitivities of CGs [ppm] c sen_NV sensitivities of NVs [ppm] c c----------------------------------------------------------------------- c LOG: c----------------------------------------------------------------------- c c 03/22/06 --bkoo-- Original development c 06/20/06 --bkoo-- 2-product model for isoprene (Henze & Seinfeld, 2006) c 08/23/13 --bkoo-- Added code for DDM-PM c 03/18/14 --bkoo-- Revised for benzene SOA c 08/25/16 --bkoo-- Updated for new SOAP c 04/22/19 --gy---- Updated hi/low NOx rate constants c c----------------------------------------------------------------------- c implicit none include 'soap.inc' integer :: iout real :: tempk,press,dthr real :: avgox(NOXID),cprec(NPREC_S),cCG(NSOAP),cNV(NNVOL) real :: yfac(3,NPREC_S) integer, parameter :: nrxn = 23 integer, parameter :: idO = 1, idOH = 2, idO3 = 3, & idNO3 = 4, idNO = 5, idHO2 = 6 integer, parameter :: idBNZA = 1, idTOLA = 2, & idXYLA = 3, idIVOA = 4, & idISP = 5, idTRP = 6, & idSQT = 7 real, parameter :: eps = 1.e-12 real*8, parameter :: rkpar(2,nrxn) = reshape( ! rate constant parameters c A E & (/ 2.30d-12, 190.0d0, ! rxn 1 BNZA + OH (CB6) & 1.80d-12, -340.0d0, ! rxn 2 TOLA + OH (CB6) & 1.85d-11, 0.0d0, ! rxn 3 XYLA + OH (CB6) & 2.70d-12, -360.0d0, ! rxn 4 ARO-O2 + NO (MCM v3.2 toluene) & 2.39d-13,-1300.0d0, ! rxn 5 ARO-O2 + HO2 (MCM v3.2 toluene) c & 2.70d-11, -390.0d0, ! rxn 6 ISP + OH (CB6) & 1.03d-14, 1995.0d0, ! rxn 7 ISP + O3 (CB6) & 3.03d-12, 448.0d0, ! rxn 8 ISP + NO3 (CB6) & 2.39d-12, -365.0d0, ! rxn 9 ISP-O2 + NO (CB6) & 7.43d-13, -700.0d0, ! rxn 10 ISP-O2 + HO2 (CB6) c & 1.50d-11, -449.0d0, ! rxn 11 TRP + OH (CB6) & 1.20d-15, 821.0d0, ! rxn 12 TRP + O3 (CB6) & 3.70d-12, -175.0d0, ! rxn 13 TRP + NO3 (CB6) & 2.70d-12, -360.0d0, ! rxn 14 TRP-O2 + NO (MCM v3.2 a-pinene) & 2.66d-13,-1300.0d0, ! rxn 15 TRP-O2 + HO2 (MCM v3.2 a-pinene) c & 1.97d-10, 0.0d0, ! rxn 16 SQT + OH (MCM v3.2 b-caryophyllene) & 1.20d-14, 0.0d0, ! rxn 17 SQT + O3 (MCM v3.2 b-caryophyllene) & 1.90d-11, 0.0d0, ! rxn 18 SQT + NO3 (MCM v3.2 b-caryophyllene) & 2.70d-12, -360.0d0, ! rxn 19 SQT-O2 + NO (assumed same as TRP-O2 + NO) & 2.66d-13,-1300.0d0, ! rxn 20 SQT-O2 + HO2 (assumed same as TRP-O2 + HO2) c & 1.34d-11, 0.0d0, ! rxn 21 IVOC + OH (Hodzic et al., 2016, ACP) & 2.70d-12, -360.0d0, ! rxn 22 IVOC-O2 + NO (assumed same as ARO-O2 + NO) & 2.39d-13,-1300.0d0 /), ! rxn 23 IVOC-O2 + HO2 (assumed same as ARO-O2 + HO2) & (/ 2, nrxn /) ) c c======================== DDM Begin ======================= c integer :: nfam,ifam real :: avgsen(nfam,NOXID) real :: sen_prec(nfam,NPREC) real :: sen_CG(nfam,NCG) real :: sen_NV(nfam,NNV) real :: srh(nfam) real :: ssumk,srk1(NOXID),sdprec logical :: lddm c c======================== DDM End ======================= c c [1/ppm-hr] = [cm3/molecule-sec]*3600*(6.022e23*1.e-12/8.314)*(P[Pa]/T[K]) real*8, parameter :: cf0 = 2.608d14 ! 3600*(6.022e23*1.e-12/8.314) real*8 :: cfac real :: rktmp(nrxn), rk1(NOXID), rh, rl, sumk, dprec, cfrac real :: dfac integer :: i, j, k logical, save :: lfirsttime = .true. c c-----Check consistency c ! if (lfirsttime) then ! lfirsttime = .false. ! endif c c-----Conversion factor for rate constants c cfac = cf0 * DBLE(press*100./tempk) ! [1/ppm-hr] = cfac * [cm3/molecule-sec] c c-----Rate constants in [1/ppm-hr] c do i = 1, nrxn rktmp(i) = SNGL( rkpar(1,i)*dexp(-rkpar(2,i)/DBLE(tempk))*cfac ) enddo c c-----Re-set NV product concentrations to zero c cNV = 0.0 c c-----NOx-dependent pathways for aromatics c rh = rktmp(4) * avgox(idNO ) ! ARO-O2 + NO rl = rktmp(5) * avgox(idHO2) ! ARO-O2 + HO2 rh = rh / (rh + rl) rl = 1.0 - rh yfac(1,idBNZA) = y_h(1,idBNZA)*rh + y_l(1,idBNZA)*rl yfac(2,idBNZA) = y_h(2,idBNZA)*rh + y_l(2,idBNZA)*rl yfac(3,idBNZA) = y_h(3,idBNZA)*rh + y_l(3,idBNZA)*rl yfac(1,idTOLA) = y_h(1,idTOLA)*rh + y_l(1,idTOLA)*rl yfac(2,idTOLA) = y_h(2,idTOLA)*rh + y_l(2,idTOLA)*rl yfac(3,idTOLA) = y_h(3,idTOLA)*rh + y_l(3,idTOLA)*rl yfac(1,idXYLA) = y_h(1,idXYLA)*rh + y_l(1,idXYLA)*rl yfac(2,idXYLA) = y_h(2,idXYLA)*rh + y_l(2,idXYLA)*rl yfac(3,idXYLA) = y_h(3,idXYLA)*rh + y_l(3,idXYLA)*rl c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam sen_NV(ifam,:) = 0.0 ! <- cNV reset srh(ifam) = ( rl * rktmp(4) * avgsen(ifam,idNO) & - rh * rktmp(5) * avgsen(ifam,idHO2) ) & / ( rktmp(4) * avgox(idNO) & + rktmp(5) * avgox(idHO2) ) c (rl)' = -(rh)' enddo endif c c======================== DDM End ======================= c rk1(idOH ) = rktmp(1) * avgox(idOH ) ! BNZA + OH cfrac = 1. - exp(-rk1(idOH )*dthr) dprec = cprec(idBNZA) * cfrac if (dprec.gt.eps) then cprec(idBNZA) = cprec(idBNZA) - dprec cCG(1) = cCG(1) + dprec * yfac(1,idBNZA) cCG(2) = cCG(2) + dprec * yfac(2,idBNZA) cNV(1) = cNV(1) + dprec * yfac(3,idBNZA) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srk1(idOH) = rktmp(1) * avgsen(ifam,idOH) sdprec = sen_prec(ifam,idBNZA) * cfrac & + cprec(idBNZA) * dthr * srk1(idOH) ! because cprec has been updated sen_prec(ifam,idBNZA) = sen_prec(ifam,idBNZA) - sdprec sen_CG(ifam,1) = sen_CG(ifam,1) + sdprec * yfac(1,idBNZA) & + dprec * & ( y_h(1,idBNZA)*srh(ifam) - y_l(1,idBNZA)*srh(ifam) ) sen_CG(ifam,2) = sen_CG(ifam,2) + sdprec * yfac(2,idBNZA) & + dprec * & ( y_h(2,idBNZA)*srh(ifam) - y_l(2,idBNZA)*srh(ifam) ) sen_NV(ifam,1) = sen_NV(ifam,1) + sdprec * yfac(3,idBNZA) & + dprec * & ( y_h(3,idBNZA)*srh(ifam) - y_l(3,idBNZA)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif rk1(idOH ) = rktmp(2) * avgox(idOH ) ! TOLA + OH cfrac = 1. - exp(-rk1(idOH )*dthr) dprec = cprec(idTOLA) * cfrac if (dprec.gt.eps) then cprec(idTOLA) = cprec(idTOLA) - dprec cCG(1) = cCG(1) + dprec * yfac(1,idTOLA) cCG(2) = cCG(2) + dprec * yfac(2,idTOLA) cNV(1) = cNV(1) + dprec * yfac(3,idTOLA) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srk1(idOH) = rktmp(2) * avgsen(ifam,idOH) sdprec = sen_prec(ifam,idTOLA) * cfrac & + cprec(idTOLA) * dthr * srk1(idOH) ! because cprec has been updated sen_prec(ifam,idTOLA) = sen_prec(ifam,idTOLA) - sdprec sen_CG(ifam,1) = sen_CG(ifam,1) + sdprec * yfac(1,idTOLA) & + dprec * & ( y_h(1,idTOLA)*srh(ifam) - y_l(1,idTOLA)*srh(ifam) ) sen_CG(ifam,2) = sen_CG(ifam,2) + sdprec * yfac(2,idTOLA) & + dprec * & ( y_h(2,idTOLA)*srh(ifam) - y_l(2,idTOLA)*srh(ifam) ) sen_NV(ifam,1) = sen_NV(ifam,1) + sdprec * yfac(3,idTOLA) & + dprec * & ( y_h(3,idTOLA)*srh(ifam) - y_l(3,idTOLA)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif rk1(idOH ) = rktmp(3) * avgox(idOH ) ! XYLA + OH cfrac = 1. - exp(-rk1(idOH )*dthr) dprec = cprec(idXYLA) * cfrac if (dprec.gt.eps) then cprec(idXYLA) = cprec(idXYLA) - dprec cCG(1) = cCG(1) + dprec * yfac(1,idXYLA) cCG(2) = cCG(2) + dprec * yfac(2,idXYLA) cNV(1) = cNV(1) + dprec * yfac(3,idXYLA) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srk1(idOH) = rktmp(3) * avgsen(ifam,idOH) sdprec = sen_prec(ifam,idXYLA) * cfrac & + cprec(idXYLA) * dthr * srk1(idOH) ! because cprec has been updated sen_prec(ifam,idXYLA) = sen_prec(ifam,idXYLA) - sdprec sen_CG(ifam,1) = sen_CG(ifam,1) + sdprec * yfac(1,idXYLA) & + dprec * & ( y_h(1,idXYLA)*srh(ifam) - y_l(1,idXYLA)*srh(ifam) ) sen_CG(ifam,2) = sen_CG(ifam,2) + sdprec * yfac(2,idXYLA) & + dprec * & ( y_h(2,idXYLA)*srh(ifam) - y_l(2,idXYLA)*srh(ifam) ) sen_NV(ifam,1) = sen_NV(ifam,1) + sdprec * yfac(3,idXYLA) & + dprec * & ( y_h(3,idXYLA)*srh(ifam) - y_l(3,idXYLA)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif c c-----NOx-dependent pathways for isoprene c rh = rktmp( 9) * avgox(idNO ) ! ISP-O2 + NO rl = rktmp(10) * avgox(idHO2) ! ISP-O2 + HO2 rh = rh / (rh + rl) rl = 1.0 - rh yfac(1,idISP) = y_h(1,idISP)*rh + y_l(1,idISP)*rl yfac(2,idISP) = y_h(2,idISP)*rh + y_l(2,idISP)*rl yfac(3,idISP) = y_h(3,idISP)*rh + y_l(3,idISP)*rl rk1(idOH ) = rktmp(6) * avgox(idOH ) ! ISP + OH rk1(idO3 ) = rktmp(7) * avgox(idO3 ) ! ISP + O3 rk1(idNO3) = rktmp(8) * avgox(idNO3) ! ISP + NO3 sumk = rk1(idOH ) + rk1(idO3 ) + rk1(idNO3) cfrac = 1. - exp(-sumk*dthr) dprec = cprec(idISP) * cfrac if (dprec.gt.eps) then cprec(idISP) = cprec(idISP) - dprec cCG(3) = cCG(3) + dprec * yfac(1,idISP) cCG(4) = cCG(4) + dprec * yfac(2,idISP) cNV(2) = cNV(2) + dprec * yfac(3,idISP) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srh(ifam) = ( rl * rktmp( 9) * avgsen(ifam,idNO) & - rh * rktmp(10) * avgsen(ifam,idHO2) ) & / ( rktmp( 9) * avgox(idNO) & + rktmp(10) * avgox(idHO2) ) ssumk = rktmp(6) * avgsen(ifam,idOH) & + rktmp(7) * avgsen(ifam,idO3) & + rktmp(8) * avgsen(ifam,idNO3) sdprec = sen_prec(ifam,idISP) * cfrac & + cprec(idISP) * dthr * ssumk ! because cprec has been updated sen_prec(ifam,idISP) = sen_prec(ifam,idISP) - sdprec sen_CG(ifam,3) = sen_CG(ifam,3) + sdprec * yfac(1,idISP) & + dprec * & ( y_h(1,idISP)*srh(ifam) - y_l(1,idISP)*srh(ifam) ) sen_CG(ifam,4) = sen_CG(ifam,4) + sdprec * yfac(2,idISP) & + dprec * & ( y_h(2,idISP)*srh(ifam) - y_l(2,idISP)*srh(ifam) ) sen_NV(ifam,2) = sen_NV(ifam,2) + sdprec * yfac(3,idISP) & + dprec * & ( y_h(3,idISP)*srh(ifam) - y_l(3,idISP)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif c c-----NOx-dependent pathways for monoterpenes c rh = rktmp(14) * avgox(idNO ) ! TRP-O2 + NO rl = rktmp(15) * avgox(idHO2) ! TRP-O2 + HO2 rh = rh / (rh + rl) rl = 1.0 - rh yfac(1,idTRP) = y_h(1,idTRP)*rh + y_l(1,idTRP)*rl yfac(2,idTRP) = y_h(2,idTRP)*rh + y_l(2,idTRP)*rl yfac(3,idTRP) = y_h(3,idTRP)*rh + y_l(3,idTRP)*rl rk1(idOH ) = rktmp(11) * avgox(idOH ) ! TRP + OH rk1(idO3 ) = rktmp(12) * avgox(idO3 ) ! TRP + O3 rk1(idNO3) = rktmp(13) * avgox(idNO3) ! TRP + NO3 sumk = rk1(idOH ) + rk1(idO3 ) + rk1(idNO3) cfrac = 1. - exp(-sumk*dthr) dprec = cprec(idTRP) * cfrac if (dprec.gt.eps) then cprec(idTRP) = cprec(idTRP) - dprec cCG(3) = cCG(3) + dprec * yfac(1,idTRP) cCG(4) = cCG(4) + dprec * yfac(2,idTRP) cNV(2) = cNV(2) + dprec * yfac(3,idTRP) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srh(ifam) = ( rl * rktmp(14) * avgsen(ifam,idNO) & - rh * rktmp(15) * avgsen(ifam,idHO2) ) & / ( rktmp(14) * avgox(idNO) & + rktmp(15) * avgox(idHO2) ) ssumk = rktmp(11) * avgsen(ifam,idOH) & + rktmp(12) * avgsen(ifam,idO3) & + rktmp(13) * avgsen(ifam,idNO3) sdprec = sen_prec(ifam,idTRP) * cfrac & + cprec(idTRP) * dthr * ssumk ! because cprec has been updated sen_prec(ifam,idTRP) = sen_prec(ifam,idTRP) - sdprec sen_CG(ifam,3) = sen_CG(ifam,3) + sdprec * yfac(1,idTRP) & + dprec * & ( y_h(1,idTRP)*srh(ifam) - y_l(1,idTRP)*srh(ifam) ) sen_CG(ifam,4) = sen_CG(ifam,4) + sdprec * yfac(2,idTRP) & + dprec * & ( y_h(2,idTRP)*srh(ifam) - y_l(2,idTRP)*srh(ifam) ) sen_NV(ifam,2) = sen_NV(ifam,2) + sdprec * yfac(3,idTRP) & + dprec * & ( y_h(3,idTRP)*srh(ifam) - y_l(3,idTRP)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif c c-----NOx-dependent pathways for sesquiterpenes c rh = rktmp(19) * avgox(idNO ) ! SQT-O2 + NO rl = rktmp(20) * avgox(idHO2) ! SQT-O2 + HO2 rh = rh / (rh + rl) rl = 1.0 - rh yfac(1,idSQT) = y_h(1,idSQT)*rh + y_l(1,idSQT)*rl yfac(2,idSQT) = y_h(2,idSQT)*rh + y_l(2,idSQT)*rl yfac(3,idSQT) = y_h(3,idSQT)*rh + y_l(3,idSQT)*rl rk1(idOH ) = rktmp(16) * avgox(idOH ) ! SQT + OH rk1(idO3 ) = rktmp(17) * avgox(idO3 ) ! SQT + O3 rk1(idNO3) = rktmp(18) * avgox(idNO3) ! SQT + NO3 sumk = rk1(idOH ) + rk1(idO3 ) + rk1(idNO3) cfrac = 1. - exp(-sumk*dthr) dprec = cprec(idSQT) * cfrac if (dprec.gt.eps) then cprec(idSQT) = cprec(idSQT) - dprec cCG(3) = cCG(3) + dprec * yfac(1,idSQT) cCG(4) = cCG(4) + dprec * yfac(2,idSQT) cNV(2) = cNV(2) + dprec * yfac(3,idSQT) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srh(ifam) = ( rl * rktmp(19) * avgsen(ifam,idNO) & - rh * rktmp(20) * avgsen(ifam,idHO2) ) & / ( rktmp(19) * avgox(idNO) & + rktmp(20) * avgox(idHO2) ) ssumk = rktmp(16) * avgsen(ifam,idOH) & + rktmp(17) * avgsen(ifam,idO3) & + rktmp(18) * avgsen(ifam,idNO3) sdprec = sen_prec(ifam,idSQT) * cfrac & + cprec(idSQT) * dthr * ssumk ! because cprec has been updated sen_prec(ifam,idSQT) = sen_prec(ifam,idSQT) - sdprec sen_CG(ifam,3) = sen_CG(ifam,3) + sdprec * yfac(1,idSQT) & + dprec * & ( y_h(1,idSQT)*srh(ifam) - y_l(1,idSQT)*srh(ifam) ) sen_CG(ifam,4) = sen_CG(ifam,4) + sdprec * yfac(2,idSQT) & + dprec * & ( y_h(2,idSQT)*srh(ifam) - y_l(2,idSQT)*srh(ifam) ) sen_NV(ifam,2) = sen_NV(ifam,2) + sdprec * yfac(3,idSQT) & + dprec * & ( y_h(3,idSQT)*srh(ifam) - y_l(3,idSQT)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif c c-----NOx-dependent pathways for IVOC (IVOA only) c rh = rktmp(22) * avgox(idNO ) ! IVOC-O2 + NO rl = rktmp(23) * avgox(idHO2) ! IVOC-O2 + HO2 rh = rh / (rh + rl) rl = 1.0 - rh yfac(1,idIVOA) = y_h(1,idIVOA)*rh + y_l(1,idIVOA)*rl yfac(2,idIVOA) = y_h(2,idIVOA)*rh + y_l(2,idIVOA)*rl yfac(3,idIVOA) = y_h(3,idIVOA)*rh + y_l(3,idIVOA)*rl c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srh(ifam) = ( rl * rktmp(22) * avgsen(ifam,idNO) & - rh * rktmp(23) * avgsen(ifam,idHO2) ) & / ( rktmp(22) * avgox(idNO) & + rktmp(23) * avgox(idHO2) ) enddo endif c c======================== DDM End ======================= c rk1(idOH ) = rktmp(21) * avgox(idOH ) ! IVOC + OH cfrac = 1. - exp(-rk1(idOH )*dthr) dprec = cprec(idIVOA) * cfrac if (dprec.gt.eps) then cprec(idIVOA) = cprec(idIVOA) - dprec cCG(1) = cCG(1) + dprec * yfac(1,idIVOA) cCG(2) = cCG(2) + dprec * yfac(2,idIVOA) cNV(1) = cNV(1) + dprec * yfac(3,idIVOA) c c======================== DDM Begin ======================= c if ( lddm ) then do ifam = 1, nfam srk1(idOH) = rktmp(21) * avgsen(ifam,idOH) sdprec = sen_prec(ifam,idIVOA) * cfrac & + cprec(idIVOA) * dthr * srk1(idOH) ! because cprec has been updated sen_prec(ifam,idIVOA) = sen_prec(ifam,idIVOA) - sdprec sen_CG(ifam,1) = sen_CG(ifam,1) + sdprec * yfac(1,idIVOA) & + dprec * & ( y_h(1,idIVOA)*srh(ifam) - y_l(1,idIVOA)*srh(ifam) ) sen_CG(ifam,2) = sen_CG(ifam,2) + sdprec * yfac(2,idIVOA) & + dprec * & ( y_h(2,idIVOA)*srh(ifam) - y_l(2,idIVOA)*srh(ifam) ) sen_NV(ifam,1) = sen_NV(ifam,1) + sdprec * yfac(3,idIVOA) & + dprec * & ( y_h(3,idIVOA)*srh(ifam) - y_l(3,idIVOA)*srh(ifam) ) enddo endif c c======================== DDM End ======================= c endif c c-----Save the weighting factors for NOx-dependent yields for PSAT c We assume same rates for ARO-O2 and IVOC-O2 c If not, PSAT needs to be re-worked c yfac(1,1) = rh yfac(2,1) = rl return end