*****************************************************************
*   MASTER MECHANISM V.3.0 ROUTINE NAME    -    RCOO2           *
*                                                               *
* UPDATES:                                                      *
* April 2015, Renee McVay, Caltech:                             *
*    Added possibility of forcing RCOO2 reactions in rcoo2.dat  *
* July 25 2013, Julia Lee-Taylor, NCAR:                         *
* 3rd channel added to RO2+HO2 [Orlando 2013, Hasson 2013]      *
*                                                               *
*            -- OLD COMMENT - NEED UPDATING --                  * 
*            -- OLD COMMENT - NEED UPDATING --                  * 
*            -- OLD COMMENT - NEED UPDATING --                  * 
*                                                               *
*   INCLUDE      general.h includes all information about glo-  *
*                          variables: CUT_OFF, ALFA, and func-  *
*                          tional groups                        *
*                                                               *   
*   COMMON BLOCK  CHMLST - Information about the chemicals in   *
*                          the dictionary                       *
*                 CHMDAT - Information about the groups and the *
*                          bond-matrix of the chemical in RDCT  *
*                                                               *
*    LOCAL CONSTANTS...                                         *
*                                                               *
*                                                               *
*    LOCAL VARIABLES...                                         *
*                                                               *
*    INTERNAL:                                                  *
*                                                               *
*    TGROUP     contents the groups of RDCT                     *
*    TBOND      contents the bond-matrix of RDCT                *
*    PCHEM      temporary variable names for products           *
*    FLAG,NR,NP flag indicates how many products (NR) exit      *
*    IND        index of chemical in CHMLST - is 0, if not ex.  *
*    NDB        number of double bonds in chain                 *
*    NCH, ICH   number and index of non-identical products      *
*    I,J,K      DO-LOOPs indices                                *
*    RATE       reaction rates of product channels              *
*               group:               2.5e-15 molecule per sec   *
*    TOTR       total sum of rates of reactions with NO3 & RDCT *
*    RTOT       total rate for reaction of double-bonded carbon *
*    FRACT      fraction deduced by position of d-bonded carbon *
*    PP,SS      products and stoichiometric coefficients after  *
*               fragmentation of excited Criegee radicals       *
*                                                               *
*    OUTPUT:                                                    *
*                                                               *
*    A1 & A4    information about type and channel of reaction  *
*    R(1-2)     reagent: input chemical and NO3                 *
*    S(N)       stoichiometric coefficients of products         *
*    P(N)       array of products: acyl radicals, HNO3 ,and co- *
*               products                                        *
*    TA         activation energy is 0.                         *
*                                                               *
* LITERATURE REFERENCES                                         *
*
*    Atkinson, R., Gas-phase tropospheric chemistry of organic
* compounds, J. Phys. Ref. Data, monograph 2, 1994.
*    Canosa-Mas, C.E., M.D. King, R. Lopez, C.J. Percival, R.P. Wayne,
* D.E. Shallcross, J.A. Pyle, and V. Daële, Is the reaction between
* CH3C(O)O2 and NO3 important in the night-time troposphere?, J.
* Chem. Soc. Faraday Trans., 92, 2211-2222, DOI: C10.1039/FT9969202211,1996.
*    Hasson, A. S.; G. S. Tyndall, J. J. Orlando, S. Singh, S. Q.
* Hernandez, S. Campbell, and Y. Ibarra, Branching ratios for the reaction
* of selected carbonyl-containing peroxy radicals with hydroperoxy
* radicals, J. Phys. Chem. A, 116, 6264-6281, doi: 10.1021/jp211799c, 2012
*    Lesclaux, R. Combination of Peroxyl Radicals in the Gas Phase. In
* Peroxyl Radicals; Alfassi, Z. B., Ed.; Wiley:  New York, 1997; pp 81−112.
*    Orlando, J.J., and G.S. Tyndall, Laboratory studies of organic
* peroxy radical chemistry: an overview with emphasis on recent issues of
* atmospheric significance, Chem. Soc. Reviews, 41, 6294-6317, 
* doi: 10.1039/C2CS35166H, 2012
*    Tomas, A., Villenave, E., and Lesclaux, R.: Reactions of the HO2
* radical with CH3CHO and CH3CO(O2) in the gas phase, J. Phys. Chem. A,
* 105, 3505–3514, https://doi.org/10.1021/jp003762p, 2001.
*    Tyndall, G. S., Cox, R. A., Granier, C., Lesclaux, R., Moortgat, G.
* K., Pilling, M. J., Ravishankara, A. R., and Wallington, T. J.,
* Atmospheric chemistry of small organic peroxy radicals, J. Geophys.
* Res., 106, 12,157– 12,182, https://doi.org/10.1029/2000JD900746, 2001.
*
*****************************************************************
      SUBROUTINE rcoo2(rdct,bond,group,nring,brch,
     &                 dbrch,dict,namlst,
     &                 cut_off,
     &                 nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &                 nfn,namfn,chemfn,
     &                 nrcoo2kr,rcoo2krct,rcoo2kprd,rcoo2arrh,rcoo2cost)

      IMPLICIT NONE
      INCLUDE 'general.h'
      INCLUDE 'organic.h'
      INCLUDE 'common.h'

* input:
      CHARACTER(LEN=lcf) rdct
      CHARACTER(LEN=lgr) group(mca)
      INTEGER         bond(mca,mca)
      INTEGER         nring
      REAL            brch
      REAL            cut_off

      CHARACTER(LEN=lfo) rcoo2krct(mkr,2),rcoo2kprd(mkr,3)
      CHARACTER(LEN=lfo) rcoo2kprd_del(mkr,3)
      INTEGER         nrcoo2kr
      REAL            rcoo2arrh(mkr,3),rcoo2cost(mkr,3)

* input/output
      CHARACTER(LEN=ldi) dict(mni)
      REAL            dbrch(mni)
      CHARACTER(LEN=lco) namlst(mni)
      INTEGER         level
      INTEGER         stabl
      CHARACTER(LEN=lst) holdvoc(mlv)
      INTEGER         nhldvoc
      CHARACTER(LEN=lst) holdrad(mra)
      INTEGER         nhldrad
      INTEGER         nfn
      CHARACTER(LEN=lco) namfn(mfn)
      CHARACTER(LEN=lfo) chemfn(mfn)

*  internal
      CHARACTER(LEN=lfo) pchem, tempkc,p_rad,p_ol
      CHARACTER(LEN=lgr) tgroup(mca), pold, pnew
      CHARACTER(LEN=lco) coprod(mca)
      INTEGER         tbond(mca,mca)
      INTEGER         ip,i,j,k,nc,nca,cnum,onum
      REAL            rmol,rrad
      REAL            brtio
      INTEGER         np,rjg(mri,2)

      CHARACTER*1     a1, a2, a3, a4
      CHARACTER(LEN=lco) p(mnp), r(3)
      REAL            s(mnp), ar1,ar2,ar3,f298,fratio
      REAL            arrh1,arrh2,arrh3,lowarrh1,lowarrh2,lowarrh3
      REAL            fc,k0M,ki,rapk,power
      REAL            ch3o2dat(3),rco3dat(3),ro2dat(8,3)
      INTEGER         idreac, nlabel
      REAL            xlabel,folow(3),fotroe(4)
      CHARACTER(LEN=lco) copchem,copchem1,copchem2
      REAL            rdtcopchem,rdtcopchem1,rdtcopchem2 
      REAL            wf
      CHARACTER(LEN=lfo) rdckprod(mca)
      CHARACTER(LEN=lco) rdcktprod(mca,mca)
      INTEGER         nip
      REAL            sc(mca),sc_del(mca),sc_temp
      CHARACTER(LEN=lfo) chem
      INTEGER         known_species,ndel
      
* data for the RO2+RO2 reactions. Following data are for the
* self reaction. First two numbers are the rate constant and
* last number is the branching ratio for the radical channel

* data for CH3O2(from Tyndall 2001)
      DATA(ch3o2dat(i),i=1,3) /9.5E-14,  -390., 0.37/

* data for RCO3 (from Tyndall 2001 for CH3COO2)
      DATA(rco3dat(i),i=1,3)  /2.5E-12,  -500., 1.0/

* data for RO2: (from Lesclaux 97)
* itype 1 : linear primary RO2
* itype 2 : branched primary RO2
* itype 3 : alpha or beta O substitued primary RO2
* itype 4 : CH3CH(OO.)CH3
* itype 5 : secondary RO2 (C>3)
* itype 6 : alpha or beta O substitued secondary RO2
* itype 7 : Tertiary RO2
* itype 8 : alpha or beta O substitued tertiary RO2
      DATA(ro2dat(1,i),i=1,3) /5.6E-14,  -500., 0.6/
      DATA(ro2dat(2,i),i=1,3) /7.8E-15, -1500., 0.5/
      DATA(ro2dat(3,i),i=1,3) /7.1E-14, -1200., 0.6/
      DATA(ro2dat(4,i),i=1,3) /1.7E-12,  2200., 0.6/
      DATA(ro2dat(5,i),i=1,3) /1.5E-13,   370., 0.3/
      DATA(ro2dat(6,i),i=1,3) /8.4E-15, -1300., 0.3/
      DATA(ro2dat(7,i),i=1,3) /4.1E-11,  4200., 1.0/
      DATA(ro2dat(8,i),i=1,3) /3.0E-13,  1220., 1.0/

      CHARACTER(LEN=lsb) :: progname='*rcoo2*   '
      CHARACTER(LEN=ler) :: mesg

* ----------
* INITIALIZE
* ----------

      IF (wtflag.NE.0) write(6,*) '*rcoo2*'

* check if species is allowed in this routine
      IF (INDEX(rdct(lco+1:lcf),acyl_peroxy).EQ.0) THEN
        WRITE(6,'(a)') '--error--, in subroutine: roo2'
        WRITE(6,'(a)') 'this routine was called with'
        WRITE(6,'(a)') 'a species with no RCO3 function:'
        WRITE(6,'(a)') rdct(lco+1:lcf)
        STOP
      ENDIF

      copchem=' '
      copchem1=' '
      copchem2=' '
      rdtcopchem=0.
      rdtcopchem1=0.
      rdtcopchem2=0. 
      p_rad=' '
      p_ol=' '
      ip = 0

* IF RINGS EXIST remove ring-join characters from groups,
      IF (nring.gt.0) THEN
        CALL rjgrm(nring,group,rjg)
      ENDIF

* initialize
      pchem=' '
      tgroup(:) = group(:)
      coprod(:) = ' '
      p(:)=' '
      tbond(:,:) = bond(:,:)
      known_species = 0

**********************************************************************
* check if the reaction is known in the database
      i = lco + INDEX(rdct(lco+1:lcf),' ')
      chem=rdct(lco+1:i)
      nc = INDEX(chem,' ') - 1
      IF (nc.LT.1) RETURN

      DO i=1,nrcoo2kr
        IF (chem.EQ.rcoo2krct(i,1)) THEN
          known_species = 1
          WRITE(6,*) 'REACTION DE RCOO2 CONNUE'
          ndel = 0
          sc_del(:)=0
          sc_temp = 0

* WRITE REACTION
          CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

* reactant
          a1   = rdct(1:1)
          r(1) = rdct(1:lco)

          ar1 = rcoo2arrh(i,1)
          ar2 = rcoo2arrh(i,2)
          ar3 = rcoo2arrh(i,3)
          f298 = ar1*(298.**ar2)*exp(-ar3/298.) 
* second reactant
          IF (rcoo2krct(i,2).NE.' ') THEN
            r(2) = rcoo2krct(i,2)(1:6)
          ENDIF

! send the products in radchk (only carbon/radicals species)
          DO j=1,3
            IF (rcoo2kprd(i,j)(1:1).EQ.' ') EXIT
            IF (INDEX(rcoo2kprd(i,j),'.').EQ.0) CYCLE 
            IF (INDEX(rcoo2kprd(i,j),'C').EQ.0) CYCLE 

            CALL radchk(rcoo2kprd(i,j),rdckprod,rdcktprod,nip,sc)
            rcoo2kprd(i,j)=rdckprod(1)           ! keep the first product
            CALL stdchm(rcoo2kprd(i,j))
            IF (nip.EQ.2) THEN                 ! if there is a delocalisation product then ...
              ndel = ndel + 1                  ! add one product
              sc_temp = rcoo2cost(i,j)           ! store the initial stochiometric coeff of the product
              rcoo2cost(i,j)= sc_temp * sc(1)    ! adjust the coefficient with the fraction of each delocalisation product
              sc_del(ndel) = sc_temp * sc(2)   ! same
              rcoo2kprd_del(i,ndel)=rdckprod(2)  ! store the delocalisation product
              DO k=1,mca                       ! store the coproducts 
                IF (rdcktprod(1,k)(1:1).NE.' ') THEN
                  ndel = ndel + 1
                  sc_del(ndel) = sc_temp * sc(1)
                  rcoo2kprd_del(i,ndel)=rdcktprod(1,k)
                ENDIF
                IF (rdcktprod(2,k)(1:1).NE.' ') THEN
                  ndel = ndel + 1
                  sc_del(ndel) = sc_temp * sc(2)
                  rcoo2kprd_del(i,ndel)=rdcktprod(2,k)
                ENDIF
              ENDDO
            ENDIF
          ENDDO

* first product
          s(1) = rcoo2cost(i,1)
          brtio = brch * s(1)
          CALL bratio(rcoo2kprd(i,1),brtio,p(1),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
          np=1
*second product
          IF (rcoo2kprd(i,2).NE.' ') THEN
            np=np+1
            s(np) = rcoo2cost(i,2)
            brtio = brch * s(np)
            CALL bratio(rcoo2kprd(i,2),brtio,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
          ENDIF
*third product
          IF (rcoo2kprd(i,3).NE.' ') THEN
            np=np+1
            s(np) = rcoo2cost(i,3)
            brtio = brch * s(np)
            CALL bratio(rcoo2kprd(i,3),brtio,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
          ENDIF
! write delocalisation products if any
          DO j=1,ndel
            IF (rcoo2kprd_del(i,j)(1:1).EQ.' ') EXIT
            np=np+1
            s(np) = sc_del(j)
            brtio = brch * s(np)
            CALL bratio(rcoo2kprd_del(i,j),brtio,p(np),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
          ENDDO
* write reaction
          CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &               f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

        ENDIF
      ENDDO
      IF (known_species.ne.0) RETURN

**********************************************************************

* locate peroxy_acyl group and count number of carbons:
      nc = INDEX(rdct(lco+1:lcf),' ') - 1
      nca = cnum(rdct(lco+1:lcf),nc)+onum(rdct(lco+1:lcf),nc)
      DO i=1,mca
         IF (INDEX(group(i),acyl_peroxy).NE.0) ip = i
      ENDDO
* -------------------
* reaction with HO2
* -------------------
* TO KEEP ONLY HIGH-NOx ways 
      IF (high_NOxfg.GT.0)      GOTO 98

* rate constant for CH3COO2 is used for all RCO3 radicals.
* Note that data for C2H5COO2 (see the review paper from
* Lesclaux, 1996) is similar to the data for CH3CO3. 
* Values are from Tomas, Villenave and Lesclaux, 2001
* J Phys Chem V105 p 3505
      arrh1 = 6.4E-13 
      arrh2 = 0.
      arrh3 = -925

      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

      ar1 = arrh1 
      ar2 = arrh2 
      ar3 = arrh3
      f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)

* to estimate the maximum production of RCO(OOH) and RCO(OH) 
* (dbrch), we assumed :
* [NO]min = 10ppt = 2.46E8 molec.cm-3
* [HO2]max = 1E9 molec.cm-3
* maximum is at 298K 1atm

* 3 channels based on:
* Orlando, J.J., and G.S. Tyndall, 2012: Chem. Soc. Reviews, 41, 6294-6317,
* doi: 10.1039/C2CS35166H.
* Hasson, A. S.; G. S. Tyndall, J. J. Orlando, S. Singh, S. Q. Hernandez,
* S. Campbell, and Y. Ibarra, 2012: J. Phys. Chem. A, 116, 6264-6281,
* doi: 10.1021/jp211799c.

* peracid channel
* ---------------

* change (OO.) to (OOH)
      pold = alkyl_peroxy
      pnew = hydro_peroxide
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild and rename:
      CALL rebond(tbond,tgroup,pchem,nring)
      CALL stdchm(pchem)

      s(1) = 0.45
      brtio = brch * s(1) / 3.

      CALL bratio(pchem,brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      
      IF (rdtcopchem.GT.0.) THEN
        WRITE(6,*)'WARNING from rcoo2.f'
        WRITE(6,*) 'Hydroperoxyde replaced by lumped species'
        WRITE(99,*) 'rcoo2',pchem !STOP
      ENDIF        

* alkoxy channel: RCO(OO.) > R(.) + OH + CO2
* -----------------------

* change (OO.) to (O)
      pold = alkyl_peroxy
      pnew = alkoxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild and rename:
      CALL rebond(tbond,tgroup,tempkc,nring)
      CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
      p_rad = rdckprod(1)
      coprod(:) = rdcktprod(1,:)
      CALL stdchm(p_rad)

      s(2) = 0.4
      brtio = brch * s(2) / 3.

      CALL bratio(p_rad,brtio,p(2),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

* OH coproduct
      s(3) = 0.4
      p(3) = 'HO    '
      s(4) = 0.4
      p(4) = 'CO2   '

* carboxylic acid channel
* -----------------------

* change (OO.) to (OH)
      pold = alkyl_peroxy
      pnew = hydroxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild and rename:
      CALL rebond(tbond,tgroup,pchem,nring)
      CALL stdchm(pchem)

      s(5) = 0.15
      brtio = brch * s(5) / 3.

      CALL bratio(pchem,brtio,p(5),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

* ozone coproduct
      s(6) = 0.15
      p(6) = 'O3   '

* -----------------------
* add coproducts from replacement by lumped species
      IF (rdtcopchem.GT.0.) THEN
        np = 4
        s(np) = rdtcopchem
        p(np) = copchem
      ENDIF

* last product : remove the counters
*      s(4) = -1.0
*      p(4) = 'XACO3'

* reactant
      a1   = rdct(1:1)
      r(1) = rdct(1:lco)
      r(2) = 'HO2  '

      fratio=1. 
      CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      IF (wtflag.NE.0) WRITE(6,*) 'done rco3 + ho2'
       
* TO KEEP ONLY ZERO-NOx ways 
      IF (zero_NOxfg.GT.0) GOTO 89

* -------------------
* reaction with NO
* -------------------
98     CONTINUE

* rate constant for CH3COO2 is used for all RCO3 radicals.
* Note that data for C2H5COO2 (see the review from
* Atkinson, 1999) is similar to the data for CH3CO3 (in
* spite of a slight increase). 
* Values are from Tyndall et al, 2001 (J. Geophysical Research)
      arrh1 = 8.1E-12
      arrh2 = 0.
      arrh3 = -270

      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

* change (OO.) to (O.)
      pold = alkyl_peroxy
      pnew = alkoxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild, check and rename:
      CALL rebond(tbond,tgroup,tempkc,nring)
      CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
      pchem = rdckprod(1)
      coprod(:) = rdcktprod(1,:)
      CALL stdchm(pchem)
      s(1) = 1.0

      CALL bratio(pchem,brch,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      IF (wtopeflag.EQ.1) THEN
        write(10,*)'********* acylper ','G',rdct(1:6),' *******'
        write(10,22)'conv NO -> NO2 /    ','G',p(1)
      ENDIF        
22    FORMAT(A21,2X,A1,A6)

* second product (for di acid)
      np = 1
      IF (nip.EQ.2) THEN
        np = np + 1
        s(1) = sc(1)
	s(np) = sc(2)	
        CALL bratio(rdckprod(2),brch,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
      ENDIF

* second product is NO2
      np = np+1
      s(np) = 1.0
      p(np) = 'NO2  '

* third product : remove the counters
*      s(3) = -1.0
*      p(3) =  'XACO3'

* other products are coproducts linked to p(1) (i.e. pchem)
c      np = 3
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(1)
	    p(np) = coprod(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (rdcktprod(2,i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = sc(2)
            p(np) = rdcktprod(2,i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (rdtcopchem.GT.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rdtcopchem
        p(np) = copchem
        IF (wtopeflag.EQ.1) 
     &        write(10,'(23X,A1,f5.3,1X,A1,A6)')'#',s(np),'G',p(np)
      ENDIF

* write:
      a1   = rdct(1:1)
      r(1) = rdct(1:lco)
      r(2) = 'NO   '

      ar1 = arrh1 
      ar2 = arrh2 
      ar3 = arrh3
      f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)
      fratio=1. 
      CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      IF (wtflag.NE.0) WRITE(6,*) 'done rcoo2 + NO'
      IF (wtopeflag.EQ.1) WRITE(10,*)'end'

* -------------------
* reaction with NO2
* -------------------

* rate constant:
* rate constant  for CH3COO2 is used for every 2 carbons species.
* Note that this reaction is in the fall-off regime. Rate constant
* are from Tyndall, 2001 (J. G.R.)
* For species with more than 2 carbons, the high pressure rate
* constant given for CH3COO2 is used, as recommended by Atkinson,
* 1994 (J. Phys. Ref. Data, monograph 2).
* For the fall off reaction, rate constant are given with the following
* expression : k=arrh1 * (T/300)**arrh2 * exp(-arrh3/T)

* low pressure rate constant
      lowarrh1 = 8.5E-29
      lowarrh2 = -6.5
      lowarrh3 = 0.

* high pressure rate constant
      arrh1 = 1.1E-11
      arrh2 = -1.0
      arrh3 = 0.

* broadening factor
      fc= 0.6

      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

* change (OO.) to (OONO2)
      pold = alkyl_peroxy
      pnew = peroxy_nitrate
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild and rename:
      CALL rebond(tbond,tgroup,pchem,nring)
      CALL stdchm(pchem)

      s(1) = 1.0
      CALL bratio(pchem,brch,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

* rien n'est prevu pour le coproduit PARC5 s'il est produit ici
      IF (rdtcopchem.GT.0.) THEN
        WRITE(6,*)'WARNING from rcoo2.f'
        WRITE(6,*) 'PAN replaced by lumped species'
        WRITE(99,*) 'rcoo2',pchem !STOP
      ENDIF        

* second product : remove the counters
*      s(2) = -1.
*      p(2) =  'XACO3'
      
* reactant:
      a1   = rdct(1:1)
      r(1) = rdct(1:lco)
      r(2) = 'NO2  '      

      IF (nca.eq.2) THEN
        r(3)='(+M)'
        s(3) = 1.
        p(3) =  '(+M)'
        ar1 = arrh1 
        ar2 = arrh2 
        ar3 = arrh3

* compute reaction rate at 298K and atmospheric pressure
        k0M=lowarrh1*((298./300.)**lowarrh2)*exp(-lowarrh3/298.)*2.45E19
        ki=arrh1*((298./300.)**arrh2)*exp(-arrh3/298.)
        rapk=k0M/ki
        power=1./(1.+log10(rapk)*log10(rapk))
        f298 = (k0m/(1.+rapk))*(fc**power)
        fratio=1. 

c        CALL rxwrit2(a1,a2,a3,a4,r,s,p,
c     &                ar1,ar2,ar3,f298,fratio,15)
c        WRITE(15,'(A6,E8.2,F5.1,F7.0,A1)')
c     &          'LOW  /',lowarrh1,lowarrh2,lowarrh3,'/'
c        WRITE(15,'(A6,F3.1,A10)')
c     &          'TROE /',fc, ' 0. 0. 0./'
c     
c        WRITE(16,'(A6,E8.2,F5.1,F7.0,A1)')
c     &          'LOW  /',lowarrh1,lowarrh2,lowarrh3,'/'
c        WRITE(16,'(A6,F3.1,A10)')
c     &          'TROE /',fc, ' 0. 0. 0./'
     
* write out 
* fall off reaction are given with idreac=3
        idreac=3
        folow(1)=lowarrh1
        folow(2)=lowarrh2
        folow(3)=lowarrh3
        fotroe(1)=fc
        CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &               f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
     
      
      ELSE
        ar1 = arrh1*((1/300.)**arrh2) 
        ar2 = arrh2 
        ar3 = arrh3
        f298 = ar1*(298.**ar2)*exp(-ar3/298.)
        fratio=1. 

        CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &               f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

      ENDIF

* -------------------
* reaction with NO3
* -------------------

* rate constant for CH3COO2 is used for all RCO3 radical.
* Value are from Canosa-Mas (J. Chem. Soc. Faraday Trans.,
* 2211-2222,1996)
      arrh1 = 5.0E-12
      arrh2 = 0.
      arrh3 = 0.

      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)

* write net NO3 reaction rate
      WRITE(71,*) rdct(1:6),' ',arrh1,arrh2,arrh3 
     &               ,' ',chem(1:index(chem,' '))

* change (OO.) to (O.)
      pold = alkyl_peroxy
      pnew = alkoxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild, check and rename:
      CALL rebond(tbond,tgroup,tempkc,nring)
      CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
      pchem = rdckprod(1)
      coprod(:) = rdcktprod(1,:)
      CALL stdchm(pchem)

      s(1) = 1.0
      brtio = brch * 0.01
      CALL bratio(pchem,brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

* second product (for di acid)
      np = 1
      IF (nip.EQ.2) THEN
        np = np + 1
        s(1) = sc(1)
        s(np) = sc(2)
        CALL bratio(rdckprod(2),brch,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
      ENDIF
* second product is NO2
      np = np + 1
      s(np) = 1.0
      p(np) = 'NO2  '

* third product : remove the counters
*      s(3) = -1.0
*      p(3) =  'XACO3'

* other products are coproducts linked to p(1) (i.e. pchem)
      np = 3
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(1)
            p(np) = coprod(i)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (rdcktprod(2,i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = sc(2)
            p(np) = rdcktprod(2,i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (rdtcopchem.GT.0.) THEN
         CALL addprod(np,progname,rdct(lco+1:lcf))
         s(np) = rdtcopchem
         p(np) = copchem
       ENDIF

* write:
      a1   = rdct(1:1)
      r(1) = rdct(1:lco)
      r(2) = 'NO3  '

      ar1 = arrh1 
      ar2 = arrh2 
      ar3 = arrh3
      f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)
      fratio=1. 

      CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      IF (wtflag.NE.0) WRITE(6,*) 'done rcoo2 + NO3'

* TO KEEP ONLY HIGH-NOx ways 
89    IF (high_NOxfg.GT.0)      GOTO 99       
* -------------------------------
* RO2+RO2 reactions
* -------------------------------
c      wf=0.001
      wf=0.02
       
* find the various products linked to the species (RC(O)O.,
* RC(O)OH)

      p_rad  = ' '
      p_ol   = ' '

      pold = alkyl_peroxy
      pnew = alkoxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)
      CALL rebond(tbond,tgroup,tempkc,nring)
      CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
      p_rad = rdckprod(1)
c      IF (nip.NE.1) STOP 'rcoo2.f 3'
      coprod(:) = rdcktprod(1,:)
      CALL stdchm(p_rad)

      pold = alkyl_peroxy
      pnew = hydroxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)
      CALL rebond(tbond,tgroup,p_ol,nring)
      CALL stdchm(p_ol)
      
* RO2 = CH3(OO.)  (NAME = CH3O2)
* ===============================

* compute the branching ratio
      rmol=(1.-ch3o2dat(3))/2.
      rrad=1.-rmol

* rate constant : Data available shows that the reaction rate
* is close to 1E-11, independant of the RO2 being considered. This
* rate constant is used, instead of the classical geometrical mean
C      arrh1 =  2.*( (rco3dat(1)*ch3o2dat(itype,1)) )**0.5
C      arrh2 =  0.
C      arrh3 =  (rco3dat(2)+ro2dat(itype,2))/2. 
      arrh1 =  1.0E-11
      arrh2 =  0.
      arrh3 =  0. 

      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      ar1 = arrh1 
      ar2 = arrh2 
      ar3 = arrh3
      f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)

* reactant
      r(1) = rdct(1:lco)
c      r(2) = 'CH3O2'
      r(2) = 'MEPERO'

* radical channel
      s(1)=rrad
c      brtio = brch * s(1) * f298 * 1E9 / 
c     &      (f298*1E9 + 7.8E-12*exp(300./298.)*2.46E8)
      brtio = brch*wf
      CALL bratio(p_rad,brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

c      s(2)=rrad
c      p(2) = 'CH3O '

* "H received" from CH3O2 channel
c      s(3)=rmol
      s(2)=rmol
c      brtio = brch * s(2) * f298 * 1E9 / 
c     &      (f298*1E9 + 7.8E-12*exp(300./298.)*2.46E8)
      brtio = brch * s(2) * wf

      CALL bratio(p_ol,brtio,p(2),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

c      s(4)=rmol
c      p(4) = 'CH2O '

* add coproducts from the radical channel (if any)
c      np=4
* second product (for di acid)
      np = 2
      IF (nip.EQ.2) THEN
        np = np + 1
c        s(1) = sc(1) * rrad
	s(np) = sc(2) * rrad
        s(1) = rrad - s(np)
        CALL bratio(rdckprod(2),brch,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
      ENDIF
c      np=2
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(1)
            p(np) = coprod(i)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (rdcktprod(2,i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = rrad * sc(2)
            p(np) = rdcktprod(2,i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (rdtcopchem1.GT.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rdtcopchem1
        p(np) = copchem1
      ENDIF
      IF (rdtcopchem2.GT.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rdtcopchem2
        p(np) = copchem2
      ENDIF

* last product : remove the counters
c      np=np+1
c      IF (np.GT.mnp) then
c         WRITE(6,'(a)') '--error--2 in rcoo2'
c         WRITE(6,'(a)') 'np is greater than mnp'
c         WRITE(6,'(a)') '(too much product in the reaction)'
c         WRITE(99,*) 'rcoo2',pchem !STOP
c      ENDIF
*      s(np) = -1.
*      p(np) =  'XACO3'

* write the reaction
      fratio=1. 

      CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      IF (wtflag.NE.0) write(6,*) 'done CH3OO'

* with other RO2 
* ================

      DO 200 j=1,8

* species for which the "self" reaction rate are lower than
* 1E-14 are here assumed to be unreactive with respect to the
* RO2+RO2 reaction. This concerns CH3CH(OO)CH3 and the tertiary
* radicals.  
        IF (j.eq.4) GOTO 200
        IF (j.eq.7) GOTO 200
        IF (j.eq.8) GOTO 200

* compute the branching ratio
        IF (ro2dat(j,3).ne.1.) THEN
          rmol=(1.-ro2dat(j,3))/2.
          rrad=1.-rmol
        ELSE
          rmol=0.
          rrad=1.
        ENDIF
        IF (wtflag.NE.0) WRITE(6,*) ro2dat(j,3), rmol,rrad

* rate constant : Data available shows that the reaction rate
* is close to 1E-11, independant of the RO2 being considered. This
* rate constant is used, instead of the classical geometrical mean
C        arrh1 =  2.*( (rco3dat(1)*ro2dat(j,1)) )**0.5
C        arrh2 =  0.
C        arrh3 =  (rco3dat(2)+ro2dat(j,2))/2. 
        arrh1 =  1.0E-11
        arrh2 =  0.
        arrh3 =  0. 

        CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &               f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
        ar1 = arrh1 
        ar2 = arrh2 
        ar3 = arrh3
        f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)

* reactant
        r(1) = rdct(1:lco)

c        r(2)='EXTRA'
c        IF (j.EQ.1) nlabel= 401
c       IF (j.EQ.2) nlabel= 402
c       IF (j.EQ.3) nlabel= 403
c        IF (j.EQ.4) nlabel= 404
c        IF (j.EQ.5) nlabel= 405
c        IF (j.EQ.6) nlabel= 406
c        IF (j.EQ.7) nlabel= 407
c        IF (j.EQ.8) nlabel= 408

c        r(2)=' '
        IF (j.EQ.1) r(2)='PERO1'
        IF (j.EQ.2) r(2)='PERO2'
        IF (j.EQ.3) r(2)='PERO3'
        IF (j.EQ.4) r(2)='PERO4'
        IF (j.EQ.5) r(2)='PERO5'
        IF (j.EQ.6) r(2)='PERO6'
        IF (j.EQ.7) r(2)='PERO7'
        IF (j.EQ.8) r(2)='PERO8'

* radical channel
        s(1)=rrad
c        brtio = brch * s(1) * f298 * 1E9 / 
c     &      (f298*1E9 + 7.8E-12*exp(300./298.)*2.46E8)
        brtio = brch * wf 
        CALL bratio(p_rad,brtio,p(1),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)

* "H received from counter" channel
        np=1
        IF (rmol.gt.0.) THEN
          CALL addprod(np,progname,rdct(lco+1:lcf))
          s(np) = rmol
c          brtio = brch * s(np) * f298 * 1E9 / 
c     &      (f298*1E9 + 7.8E-12*exp(300./298.)*2.46E8)
          brtio = brch * s(np) * wf
          CALL bratio(p_ol,brtio,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
        ENDIF

* second product (for di acid)
c      np = 2
      IF (nip.EQ.2) THEN
        np = np + 1
c        s(1) = sc(1) * rrad
	s(np) = sc(2) * rrad
        s(1) = rrad - s(np)
        CALL bratio(rdckprod(2),brch,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
      ENDIF

* add coproducts from the radical channel (if any)
        DO i=1,mca
           IF (coprod(i)(1:1).NE.' ') THEN
              CALL addprod(np,progname,rdct(lco+1:lcf))
              s(np) = s(1)
              p(np) = coprod(i)
           ENDIF
        ENDDO
        DO i=1,mca
          IF (rdcktprod(2,i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = sc(2) * rrad
            p(np) = rdcktprod(2,i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
          ENDIF
        ENDDO
        IF (rdtcopchem1.GT.0.) THEN
           CALL addprod(np,progname,rdct(lco+1:lcf))
           s(np) = rdtcopchem1
           p(np) = copchem1
         ENDIF
        IF (rdtcopchem2.GT.0.) THEN
           CALL addprod(np,progname,rdct(lco+1:lcf))
           s(np) = rdtcopchem2
           p(np) = copchem2
         ENDIF

* write the reaction
        fratio=1. 

* write out - extra reaction => idreac=2 (nlabel was set above)
c        idreac=2
        idreac=0
        nlabel=0
        CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
        IF (wtflag.NE.0) write(6,*) 'done rcoo2 reaction with counter 
     &        number',j

200   CONTINUE

* with RCO3 (acyl peroxy counters):
* =================================

* compute the branching ratio
       rrad=1.

* rate constant :  2*(self reaction rate)
      arrh1 =  2.*rco3dat(1)
      arrh2 =  0.
      arrh3 =  rco3dat(2)
        
      CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      ar1 = arrh1 
      ar2 = arrh2 
      ar3 = arrh3
      f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)

* reactant
      r(1) = rdct(1:lco)
c      r(2) = 'XACO3'
c      r(2) = 'EXTRA'
      r(2) = 'PERO9'
c      nlabel=409

* radical channel
      s(1)=rrad
c      brtio = brch * s(1) * f298 * 1E9 / 
c     &      (f298*1E9 + 7.8E-12*exp(300./298.)*2.46E8)
      brtio = brch * wf 
      CALL bratio(p_rad,brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

* second product (for di acid)
      np = 1
      IF (nip.EQ.2) THEN
        np = np + 1
c        s(1) = sc(1) * rrad
	s(np) = sc(2) * rrad
        s(1) = rrad - s(np)
        CALL bratio(rdckprod(2),brch,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
      ENDIF

* add coproducts from the radical channel (if any)
c      np=1
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(1)
            p(np) = coprod(i)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (rdcktprod(2,i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = sc(2) * rrad
            p(np) = rdcktprod(2,i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (rdtcopchem.GT.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rdtcopchem
        p(np) = copchem
      ENDIF

* write the reaction
      fratio=1. 

* write out - extra reaction => idreac=2 (nlabel was set above)
c      idreac=2
      idreac=0
      nlabel=0
      CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3,
     &             f298,fratio,idreac,nlabel,xlabel,folow,fotroe)
      IF (wtflag.NE.0) write(6,*) 'done RCO3'

* end of RO2+RO2 reaction
* end RCOO2
99    CONTINUE
      RETURN
      END

