*****************************************************************
*   MASTER MECHANISM V.3.0 ROUTINE NAME    -    RO2             *
*                                                               *
* UPDATES:                                                      *
* June 16 2015, Julia Lee-Taylor, NCAR:                         *
* RO2+HO2 2nd channel restricted to beta-carbonyls ONLY         *
* July 25 2013, Julia Lee-Taylor, NCAR:                         *
* 2nd 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 *
*    DFRACT     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.                         *
*                                                               *
*****************************************************************
      SUBROUTINE ro2(rdct,bond,group,nring,brch,temp,tbm,
     &               dbrch,dict,namlst,
     &               cut_off,
     &               nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &               nfn,namfn,chemfn,
     &              nro2kr,ro2krct,ro2kprd,ro2arrh,ro2cost)
      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),nring
      REAL            cut_off
      REAL            temp
      REAL            tbm

      CHARACTER(LEN=lfo) ro2krct(mkr,2),ro2kprd(mkr,3),
     &                   ro2kprd_del(mkr,3)
      INTEGER         nro2kr
      REAL            ro2arrh(mkr,3),ro2cost(mkr,3)

* input/output
      REAL            brch
      CHARACTER(LEN=ldi) dict(mni)
      REAL            dbrch(mni)
      CHARACTER(LEN=lco) namlst(mni)

      INTEGER          nfn
      CHARACTER(LEN=lco)  namfn(mfn)
      CHARACTER(LEN=lfo)  chemfn(mfn)

* internal
      CHARACTER(LEN=lfo) pchem(mca), tempkc, tempkcp 
      CHARACTER(LEN=lfo) p_onyl, p_ol, p_rad(2)
      CHARACTER(LEN=lgr) tgroup(mca), tempkg, pold, pnew
      CHARACTER(LEN=lco) coprod(mca),coprod2(mca),coprod_alk(mca)
      INTEGER         tbond(mca,mca)
      INTEGER         np,ip,i,j,k,l,itype,nc,nca,cnum,onum,ngrp
      INTEGER         ipo,icl,ibr,iyl,flag,xchk
      REAL            rab,rad,fract,brtio,nityield,check
      REAL            xc,rrad,rmol,rmol1,rmol2
      INTEGER         level
      INTEGER         stabl
      CHARACTER(LEN=lst) holdvoc(mlv)
      INTEGER         nhldvoc
      CHARACTER(LEN=lst) holdrad(mra)
      INTEGER         nhldrad
      CHARACTER*1     a1, a2, a3, a4
      CHARACTER(LEN=lco) p(mnp), r(3)
      CHARACTER*5     wrtnum
      REAL            s(mnp), ar1,ar2,ar3,f298,fratio
      REAL            arrh1,arrh2,arrh3
      REAL            ch3o2dat(3),rco3dat(3),ro2dat(8,3)
      INTEGER         idreac, nlabel
      REAL            xlabel,folow(3),fotroe(4)
      INTEGER         rjg(mri,2)
      CHARACTER(LEN=lco) copchem,copchem1,copchem2,copchem3
      REAL            rdtcopchem,rdtcopchem1,rdtcopchem2,rdtcopchem3 
      REAL            wf
      CHARACTER(LEN=lfo) rdckprod(mca)
      CHARACTER(LEN=lco) rdcktprod(mca,mca),prod_ol
      INTEGER         nip
      REAL            sc(mca),sc_del(mca),sc_temp
      CHARACTER(LEN=lfo) chem
      INTEGER         known_species,ndel
      INTEGER         :: track(mco,mca)
      INTEGER         :: trlen(mco)
      INTEGER         :: ntr
      INTEGER         :: ring(mca),rngflg
! number of atoms
      INTEGER         :: ic,ih,in,io,ir,is,ifl,ib,iclor
! neighbouring functions
      INTEGER         :: n_betaOH, n_alphaCO, n_betaONO2, n_betaOOH 
      INTEGER         :: n_gammaOH, n_betaCO, n_gammaONO2, n_gammaOOH
! chemmap output
      INTEGER         :: nfunc
      CHARACTER*1     :: nodetype(mca)      
      REAL            :: alifun(20),cdfun(20),arofun(20)
      REAL            :: mapfun(mca,3,20)
      INTEGER         :: funflg(mca)
      INTEGER         :: tabester(4,2)  ! 1= -O- side, 2= CO side
      INTEGER         :: nfcd,nfcr
      INTEGER         :: ierr
! abcde_map output
      INTEGER         :: nabcde(9), tabcde(9,mco,mca)
      
* 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.0E-10,  2200., 0.6/
      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='*ro2*   '
      CHARACTER(LEN=ler) :: mesg

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

      IF (wtflag.NE.0) write(*,*) '*ro2*'

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

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

* initialize
      tgroup(:) = group(:)
      tbond(:,:)=bond(:,:)
      coprod(:) = ' '
      coprod2(:) = ' '
      coprod_alk(:)= ' '
      pchem(:) = ' '
      copchem=' '
      copchem1=' '
      copchem2=' '
      copchem3=' '
      rdtcopchem=0.
      rdtcopchem1=0.
      rdtcopchem2=0.
      rdtcopchem3=0. 
      known_species=0.
      ro2kprd_del(:,:)=' '

**********************************************************************
* 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,nro2kr
        IF (chem.EQ.ro2krct(i,1)) THEN
          known_species = 1
          WRITE(6,*) 'REACTION DE RO2 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 = ro2arrh(i,1)
          ar2 = ro2arrh(i,2)
          ar3 = ro2arrh(i,3)
          f298 = ar1*(298.**ar2)*exp(-ar3/298.) 
* second reactant
          IF (ro2krct(i,2).NE.' ') THEN
            r(2) = ro2krct(i,2)(1:6)
          ENDIF

! send the products in radchk (only carbon/radicals species)
          DO j=1,3
            IF (ro2kprd(i,j)(1:1).EQ.' ') EXIT
            IF (INDEX(ro2kprd(i,j),'.').EQ.0) CYCLE 
            IF (INDEX(ro2kprd(i,j),'C').EQ.0) CYCLE 
            CALL radchk(ro2kprd(i,j),rdckprod,rdcktprod,nip,sc)
            ro2kprd(i,j)=rdckprod(1)           ! keep the first product
            CALL stdchm(ro2kprd(i,j))
            IF (nip.EQ.2) THEN                 ! if there is a delocalisation product then ...
              ndel = ndel + 1                  ! add one product
              sc_temp = ro2cost(i,j)           ! store the initial stochiometric coeff of the product
              ro2cost(i,j)= sc_temp * sc(1)    ! adjust the coefficient with the fraction of each delocalisation product
              sc_del(ndel) = sc_temp * sc(2)   ! same
              ro2kprd_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)
                  ro2kprd_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)
                  ro2kprd_del(i,ndel)=rdcktprod(2,k)
                ENDIF
              ENDDO
            ENDIF
          ENDDO

* first product
          s(1) = ro2cost(i,1)
          brtio = brch * s(1)
          CALL bratio(ro2kprd(i,1),brtio,p(1),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
          np=1
*second product
          IF (ro2kprd(i,2).NE.' ') THEN
            np=np+1
            s(np) = ro2cost(i,2)
            brtio = brch * s(np)
            CALL bratio(ro2kprd(i,2),brtio,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
          ENDIF
*third product
          IF (ro2kprd(i,3).NE.' ') THEN
            np=np+1
            s(np) = ro2cost(i,3)
            brtio = brch * s(np)
            CALL bratio(ro2kprd(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 (ro2kprd_del(i,j)(1:1).EQ.' ') EXIT
            np=np+1
            s(np) = sc_del(j)
            brtio = brch * s(np)
            CALL bratio(ro2kprd_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 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)
! also count other atoms, needed for new ro2+ho2 method
      call number(rdct(lco+1:lcf),nc,ic,ih,in,io,ir,is,ifl,ib,iclor)
      DO i=1,mca
         IF (INDEX(group(i),alkyl_peroxy).NE.0)  ip = i
      ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ric --- AROMATIC CHEMISTRY MAGNIFY 2016
! Check for special structures from OH addition on aromatic ring
      CALL gettrack(tbond,ip,nca,ntr,track,trlen)
      ngrp=0
      DO i=1,mca
        IF (group(i).NE.' ') ngrp=ngrp+1
      ENDDO
      DO i=1,nca 
        IF ((tbond(ip,i).NE.0).AND.(INDEX(group(i),'(OH)').NE.0)) THEN
          CALL findring(ip,i,nca,tbond,rngflg,ring)
          IF (rngflg.EQ.1) THEN
          DO j=1,ntr
            IF (trlen(j).LT.6) CYCLE
            IF ((tgroup(track(j,2))(1:2).EQ.'Cd').AND.
     &      (tgroup(track(j,3))(1:2).EQ.'Cd').AND.
     &      (tgroup(track(j,4))(1:2).EQ.'Cd').AND.
     &      (tgroup(track(j,5))(1:2).EQ.'Cd').AND.
     &      (track(j,6).EQ.i)) THEN

! 1st reaction is a ring closure from with a -O--O-bridge within the carbon
! ring, followed by stabilisation
! add new groups :
              IF (wtflag.NE.0) WRITE(6,*) 'RO2 ring closure'
              arrh1 = 0.75
              arrh2 = 0.
              arrh3 = 0.

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

              tgroup(ngrp+1)='-O-'  
              tgroup(ngrp+2)='-O-'
              pold='(OO.)'
              pnew=' '
              CALL swap(group(ip),pold,tgroup(ip),pnew)
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,4)),pold,tgroup(track(j,4)),pnew)
              nc=INDEX(tgroup(track(j,4)),' ')
              tgroup(track(j,4))(nc:nc)='.'
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,5)),pold,tgroup(track(j,5)),pnew)
              CALL setbond(tbond,ngrp+1,ip,3)
              CALL setbond(tbond,ngrp+1,ngrp+2,3)
              CALL setbond(tbond,ngrp+2,track(j,5),3)
              CALL setbond(tbond,track(j,4),track(j,5),1)
            
              CALL rebond(tbond,tgroup,tempkc,nring)
              CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
              CALL stdchm(rdckprod(1))
              pchem(1)=rdckprod(1)
              s(1)=1
              IF (nip.GT.1) THEN
                CALL stdchm(rdckprod(2))
                pchem(2)=rdckprod(2)
                s(2) = sc(2)
                s(1) = sc(1)
              ENDIF
              
              brtio = brch  ! force to highest yield
              CALL bratio(pchem(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
              np = 1
              IF (nip.EQ.2) THEN
              np = np + 1
              CALL bratio(pchem(np),brtio,p(np),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
              ENDIF

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

              ar1 = arrh1
              ar2 = arrh2
              ar3 = arrh3
              f298 = arrh1*(298.**arrh2)*exp(-arrh3/298.)
              fratio=1.
!              WRITE(6,*) 'p(1)',p(1),s(1)
!              WRITE(6,*) 'p(2)',p(2),s(2)

              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 ro2 ring closure'

              tbond(:,:)=bond(:,:)
              tgroup(:)=group(:)

! OO ring closure followed by O--O breaking, leading to an epoxy
! function and an alkoxy radical.
              IF (wtflag.NE.0) WRITE(6,*) 'post aromatic chemistry'
!              arrh1 = 0.125
!              arrh2 = 0.
!              arrh3 = 0.

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

              tgroup(ngrp+1)='-O-'
              pold='(OO.)'
              pnew='(O.)'
              CALL swap(group(ip),pold,tgroup(ip),pnew)
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,4)),pold,tgroup(track(j,4)),pnew)
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,5)),pold,tgroup(track(j,5)),pnew)
              CALL setbond(tbond,track(j,4),track(j,5),1)
              CALL setbond(tbond,track(j,4),ngrp+1,3)
              CALL setbond(tbond,track(j,5),ngrp+1,1)

              CALL rebond(tbond,tgroup,tempkc,nring)
              CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
              CALL stdchm(rdckprod(1))
              pchem(1)=rdckprod(1)
!              s(1)=1

!              brtio = brch  ! force to highest yield
!              CALL bratio(pchem(1),brtio,p(1),
!     &            dbrch,dict,namlst,
!     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
!     &            nfn,namfn,chemfn)
!              np = 1

* reactant
!              a1   = rdct(1:1)
!              r(1) = rdct(1:lco)
!              r(2) = '     '
!
!              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 epoxy formation'

              tbond(:,:)=bond(:,:)
              tgroup(:)=group(:)

! OO ring closure followed byt O--O breaking, leading to an epoxy
! function and an alkoxy radical.
              IF (wtflag.NE.0) WRITE(6,*) 'post aromatic chemistry'

              tgroup(ngrp+1)='-O-'
              pold='(OO.)'
              pnew=' '
              CALL swap(group(ip),pold,tgroup(ip),pnew)
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,2)),pold,tgroup(track(j,2)),pnew)
              pold='Cd'
              pnew='C'
              CALL swap(group(track(j,5)),pold,tgroup(track(j,5)),pnew)

              nc = INDEX(tgroup(track(j,5)),' ')
              tgroup(track(j,5))(nc:nc+3)='(O.)'

              CALL setbond(tbond,track(j,2),track(j,3),1)
              CALL setbond(tbond,track(j,3),track(j,4),2)
              CALL setbond(tbond,track(j,4),track(j,5),1)
              CALL setbond(tbond,ip,ngrp+1,3)
              CALL setbond(tbond,track(j,2),ngrp+1,3)

              CALL rebond(tbond,tgroup,tempkc,nring)
              CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
              CALL stdchm(rdckprod(1))
              pchem(2)=rdckprod(1)

              DO k=1,2
                IF (pchem(2).EQ.pchem(1)) THEN
                  arrh1 = 0.25
                  IF (k.EQ.2) CYCLE
                ELSE
                  arrh1 = 0.125
                ENDIF
                arrh2 = 0.
                arrh3 = 0.

!                IF (pchem(2).EQ.pchem(1)) CYCLE

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

                brtio = brch  ! force to highest yield
                CALL bratio(pchem(k),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

                s(1)=1

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

                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 epoxy formation'
              ENDDO

              tbond(:,:)=bond(:,:)
              tgroup(:)=group(:)

              RETURN

!            WRITE(6,*) sc(1),rdckprod(1)
!            WRITE(6,*) sc(2),rdckprod(2)
!            DO l=1,nca
!              WRITE(6,*) l,group(l)
!            ENDDO
!            WRITE(6,*) (track(j,k),k=1,trlen(j))
!            WRITE(6,*) '----------------'
!            DO l=1,ngrp
!              WRITE(6,*) l,tgroup(l)
!            ENDDO
!            WRITE(6,*) '----------------'

            ENDIF
          ENDDO
!          STOP
          ENDIF

        ENDIF
      ENDDO


      tbond(:,:)=bond(:,:)
      tgroup(:)=group(:)

* identify type of RO2:
* 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

* count for RO2 class (primary, secondary, ...), branched
* structure in beta, polar subtitution in alpha or beta,
* carbonyl in alpha or beta (but a-carbonyls go to rcoo2.f)
      icl=0
      ibr=0
      ipo=0
      iyl=0
      DO i=1, mca
         IF (bond(ip,i).ne.0) THEN
            icl=icl+1
            DO j=1, mca
              IF ( (bond(i,j).ne.0).AND.(j.ne.ip) ) ibr=ibr+1
            ENDDO
            IF (INDEX(group(i),'O').NE.0) ipo=1
            IF (INDEX(group(i),'N').NE.0) ipo=1
            IF (group(i).EQ.carbonyl.OR.
     &          group(i).EQ.aldehyde.OR.
     &          group(i).EQ.ketone) iyl=1
         ENDIF
      ENDDO

* assign the type. Note that substitution at the carbon bearing
* the peroxy group is not tested above. This is done below if
* the number of H is not the number expected.
      itype=0
      IF (icl.EQ.1) THEN
         itype=1
         IF (INDEX(group(ip),'CH2').EQ.0) itype=3
         IF (ibr.gt.1) itype=2
         IF (ipo.ne.0) itype=3
      ELSE IF (icl.EQ.2) THEN
         itype = 4
         IF (INDEX(group(ip),'CH').EQ.0) itype=6
         IF (nca.gt.3) itype=5
         IF (ipo.ne.0) itype=6
      ELSE IF (icl.EQ.3) THEN
         itype = 7
         IF (ipo.ne.0) itype=8
      ENDIF
* check that a type was assigned to the RO2
      IF (itype.eq.0) THEN
         WRITE(6,'(a)') '--warning, (stop) in ro2'
         WRITE(6,'(a)') 'the following species has no type'
         WRITE(6,'(a)') '(e.g. primary, secondary, ...)'
         WRITE(6,'(a)') rdct(lco+1:lcf)
         WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF

* -------------------
* reaction with NO
* -------------------

* rate constant :
* It has been assumed in the past that the reaction rate decreases
* when the size of the RO2 increases (e.g. Jenkin et al., Atm. Env.,
* 81-104, 1997). However, value provided in the recent review of
* Atkinson (1999, 1997) as well as recent rate constant measured
* by Howard et al. (1996) does not clearly show such a trend. 
* RO2+NO reaction rate are based on Atkinson recommendation 
* (J. Phys. Chem. Ref. Data, 1997), i.e. k=2.7E-12*exp(360/T) for
* all RO2
      arrh1 =  2.7E-12
      arrh2 =  0.
      arrh3 = -360. 

* nitrate yield
* function nityield compute the nitrate yield for unsubstitued
* secondary RO2 as a function of the carbon number. For other
* RO2, SAPRC99 (Carter, 1999) parameterization is used, i.e.  
* carbon number is decreased by 1.5. If the function nityield
* is called with flag=0, then nitrate yield for the input 
* temperature and pressure, otherwise standard conditions are
* used (temp=298K,pres=2.46E19) 
      flag=0
      xchk=0

* non secondary (or substituted) RO2
      IF (INDEX(group(ip),'CH(OO.)').EQ.0)  xchk=1

* substituted RO2
      DO i=1,mca
        IF ( (INDEX(group(i),'O').ne.0).AND.(i.ne.ip) ) xchk=1
        IF (INDEX(group(i),'N').ne.0)                   xchk=1
        IF ((group(i)(1:8).EQ.'CH2(OH) ').OR.
     &      (group(i)(1:7).EQ.'CH(OH) ').OR.
     &      (group(i)(1:6).EQ.'C(OH) ')) THEN
           xchk=2
           GOTO 122
        ENDIF
      ENDDO
122   CONTINUE
      
      IF (xchk.eq.0) THEN
        xc=REAL(nca)
      ELSE
        xc=REAL(nca)-1.5
      ENDIF

      IF (nca.GE.3) THEN
         fract= nityield(xc,itype,flag,temp,tbm)
         IF (xchk.EQ.2) fract=fract/2
         rad = fract
         rab = (1.-fract)
      ELSE
         rad   = 0.
         rab   = 1.
      ENDIF

! forbid the nitrate formation if there is a functional group close to
! the peroxy radical (march 2015)
      j = 0
      DO i=1,nca
        IF (tbond(ip,i).EQ.3) THEN   ! check if there is an ether group next to the peroxy 
          j = 1
          EXIT
        ENDIF
      ENDDO
      IF (((group(ip)(1:9).NE.'CH2(OO.) ').AND.
     &     (group(ip)(1:8).NE.'CH(OO.) ').AND. 
     &     (group(ip)(1:7).NE.'C(OO.) ')).OR.
     &     (j.NE.0)) THEN
         rad   = 0.
         rab   = 1.
      ENDIF
      
* RO2+NO => RO+NO2 reaction
* -------------------------
      IF (wtflag.NE.0) WRITE(6,*) 'NO abstraction'

      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(1) = rdckprod(1)
      CALL stdchm(pchem(1))
      s(1) = 1.
      IF (nip.EQ.2) THEN
        pchem(2) = rdckprod(2)
        CALL stdchm(pchem(2))
        s(2) = sc(2)
        s(1) = sc(1)
      ENDIF

* check if product is substituted alkene
* (and decompose it if necessary)
      np = 0
      DO k=1,nip
        IF(INDEX(pchem(k),'CdH(').NE.0.OR.
     &     INDEX(pchem(k),'Cd(').NE.0) THEN
          DO i=1,mca
            IF(INDEX(rdcktprod(k,i),' ').eq.1) THEN
              CALL alkcheck(pchem(k),rdcktprod(k,i))
              tempkcp = rdcktprod(k,i)
              IF(INDEX(pchem(k),'.').NE.0) THEN
                tempkc = pchem(k)
                CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
                IF (nip.EQ.2) THEN
                  PRINT*,"error in ro2 + no, decompose subst. alkene"
                  PRINT*,"------------> too many products!"
                  STOP
                ENDIF
                pchem(k) = rdckprod(k)
              ENDIF
              CALL stdchm(pchem(k))
              EXIT
            ENDIF
          ENDDO
* second product is coproduct from alkcheck
* (code will stop if radchk produces additional coproducts)
          np = np + 1
          coprod_alk(np)=tempkcp
!          s(np) = 1.
!          p(np) = tempkcp
        ENDIF
      ENDDO

      coprod(:) = rdcktprod(1,:)
      coprod2(:) = rdcktprod(2,:)

c      s(1) = 1.

* brtio : local branching ratio
* brch : branching ratio for the total production by this way
* dbrch : sum of all the production ways for this species
c      brtio = rab * brch
      brtio = brch  ! force to highest yield
      CALL bratio(pchem(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      np = 1
      IF (nip.EQ.2) THEN
      np = np + 1
      CALL bratio(pchem(np),brtio,p(np),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      ENDIF
      IF (wtopeflag.EQ.1) THEN
      write(10,*)'********** peroxy ','G',rdct(1:6),':',itype,' *******'
      write(10,22)'conv NO -> NO2 /    ',rab,'G',p(1)
22     FORMAT(A21,f5.3,2X,A1,A6)
      ENDIF
* second product is NO2
      np = np + 1
      s(np) = 1.
      p(np) = 'NO2  '

* third product : remove the counters
*      s(3) = -1.
*      IF (itype.EQ.1) p(3) = 'XP1O2'
*      IF (itype.EQ.2) p(3) = 'XP2O2'
*      IF (itype.EQ.3) p(3) = 'XP3O2'
*      IF (itype.EQ.4) p(3) = 'XS1O2'
*      IF (itype.EQ.5) p(3) = 'XS2O2'
*      IF (itype.EQ.6) p(3) = 'XS3O2'
*      IF (itype.EQ.7) p(3) = 'XT1O2'
*      IF (itype.EQ.8) p(3) = 'XT2O2'

* other products are coproducts linked to p(1) (i.e. pchem(1))
c      np = 3
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (nip.EQ.2) THEN
      DO i=1,mca
         IF (coprod2(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod2(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      ENDIF
      DO i=1,mca
         IF (coprod_alk(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod_alk(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      IF (rdtcopchem.GT.0.) THEN
        WRITE(6,*)' voie NO->NO2 a donne du PARC5 (ro2.f)'
        WRITE(6,*)' rien de prevu ds operateur.f pour ca'
        WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF

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

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

      IF (rad.NE.0.) a4 = 'A'
      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 ro2 + no abstraction'


* RO2+NO => RONO2 reaction
* ------------------------

      IF (rad.NE.0.) THEN
        IF (wtflag.NE.0) print*,'NO addition'

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

* change (OO.) to (ONO2)
        pold = alkyl_peroxy  ! JMLT added
        pnew = nitrate
        CALL swap(group(ip),pold,tgroup(ip),pnew)

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

* check if product is substituted alkene
* (and decompose it if necessary)
        IF(INDEX(pchem(1),'CdH(').NE.0.OR.
     &     INDEX(pchem(1),'Cd(').NE.0) THEN
          np=0
          DO i=1,mca
            IF(INDEX(rdcktprod(1,i),' ').eq.1) THEN
              CALL alkcheck(pchem(1),rdcktprod(1,i))
              tempkcp = rdcktprod(1,i)
              IF(INDEX(pchem(1),'.').NE.0) THEN
                tempkc = pchem(1)
                CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
                IF (nip.EQ.2) THEN
                  PRINT*,"error in ro2+no > rono2, subst. alkene"
                  PRINT*,"------------> too many products!"
                  STOP
                ENDIF
                pchem(1) = rdckprod(1)
              ENDIF
              CALL stdchm(pchem(1))
              s(1) = 1.
              rdcktprod(1,i) = tempkcp
              EXIT
            ENDIF
* second product is coproduct from alkcheck
* (code will stop if radchk produces additional coproducts)
            np = np + 1
            coprod_alk(np)=tempkcp
!            s(np) = 1.
!            p(np) = tempkcp
          ENDDO
        ENDIF

        coprod(:) = rdcktprod(1,:)

        s(1) = 1.
        brtio = rad * brch
        CALL bratio(pchem(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
       
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (coprod_alk(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod_alk(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
       IF (wtopeflag.EQ.1) THEN
         IF (rdtcopchem.GT.0.) THEN
           s(2) = rdtcopchem
           p(2) = copchem
           write(10,23)'conv NO -> nitrate / ',rad,'G',p(1),
     &           '+',s(2)*rad,'G',p(2)
         ELSE
           write(10,22)'conv NO -> nitrate / ',rad,'G',p(1)
         ENDIF
23       FORMAT(A21,f5.3,2X,A1,A6,1X,A1,1X,f5.3,1X,A1,A6)
       ENDIF

* second product : remove the counters
*        s(2) = -1.
*        IF (itype.EQ.1) p(2) = 'XP1O2'
*        IF (itype.EQ.2) p(2) = 'XP2O2'
*        IF (itype.EQ.3) p(2) = 'XP3O2'
*        IF (itype.EQ.4) p(2) = 'XS1O2'
*        IF (itype.EQ.5) p(2) = 'XS2O2'
*        IF (itype.EQ.6) p(2) = 'XS3O2'
*        IF (itype.EQ.7) p(2) = 'XT1O2'
*        IF (itype.EQ.8) p(2) = 'XT2O2'

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

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

        a4   = 'B'
        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 ro2 + no addition'
      ENDIF
      IF (wtopeflag.EQ.1) WRITE(10,*)'end'

* -------------------
* JMLT, June 2020, RO2+NO2 sensitivity test.
* -------------------
* reaction with NO2 (adapted from NO reaction section)
* -------------------
      IF (ro2no2_fg.EQ.1) THEN

* rate constant (per John Orlando, pers. comm., June 2020):
      arrh1 =  5.0E-12
      arrh2 =  0.
      arrh3 =  0. 

* RO2+NO2 => RO2NO2 reaction
* ------------------------

      IF (rad.NE.0.) THEN
        IF (wtflag.NE.0) print*,'NO2 addition'

        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(1),nring)
        CALL stdchm(pchem(1))

* check if product is substituted alkene
* (and decompose it if necessary)
        IF(INDEX(pchem(1),'CdH(').NE.0.OR.
     &     INDEX(pchem(1),'Cd(').NE.0) THEN
          np=0
          DO i=1,mca
            IF(INDEX(rdcktprod(1,i),' ').eq.1) THEN
              CALL alkcheck(pchem(1),rdcktprod(1,i))
              tempkcp = rdcktprod(1,i)
              IF(INDEX(pchem(1),'.').NE.0) THEN
                tempkc = pchem(1)
                CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
                IF (nip.EQ.2) THEN
                  PRINT*,"error in ro2+no2 > ro2no2, subst. alkene"
                  PRINT*,"------------> too many products!"
                  STOP
                ENDIF
                pchem(1) = rdckprod(1)
              ENDIF
              CALL stdchm(pchem(1))
              s(1) = 1.
              rdcktprod(1,i) = tempkcp
              EXIT
            ENDIF
* second product is coproduct from alkcheck
* (code will stop if radchk produces additional coproducts)
            np = np + 1
            coprod_alk(np)=tempkcp
!            s(np) = 1.
!            p(np) = tempkcp
          ENDDO
        ENDIF

        coprod(:) = rdcktprod(1,:)

        s(1) = 1.
        brtio = 1.
        CALL bratio(pchem(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
       
      DO i=1,mca
         IF (coprod(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      DO i=1,mca
         IF (coprod_alk(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = 1.
            p(np) = coprod_alk(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
       IF (wtopeflag.EQ.1) THEN
         IF (rdtcopchem.GT.0.) THEN
           s(2) = rdtcopchem
           p(2) = copchem
           write(10,23)'conv NO2-> perono2 / ',rad,'G',p(1),
     &           '+',s(2)*rad,'G',p(2)
         ELSE
           write(10,22)'conv NO2-> perono2 / ',rad,'G',p(1)
         ENDIF
       ENDIF

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

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

        a4   = 'B'
        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 ro2 + no2 addition'
      ENDIF
      IF (wtopeflag.EQ.1) WRITE(10,*)'end'

      ENDIF ! ro2no2_fg
* END JMLT, June 2020, RO2+NO2 sensitivity test.

* TO KEEP ONLY HIGH-NOx ways 
      IF (high_NOxfg.GT.0) GOTO 99

* -------------------
* reaction with HO2
* -------------------
      IF (wtflag.NE.0) print*,'rxn with HO2'

      IF (ro2ho2_fg .eq. 2) then
* rate constant according to Wennberg et al., 2018
* expression is
* k = 2.82e-13*exp(1300/T)*(1-exp(-0.231*(nc + no + nn -2)))
        arrh1 = 2.82e-13*(1-exp(-0.231*(ic + io + in - 2)))
        arrh2 = 0.
        arrh3 = -1300.      
      else
* rate constant. The following expression was found in
* "advances in European Tropospheric research", Air pollution 
* research report 72 for T=298K : k=2.28E-11*(1-exp(-0.245*nca))
* where nca is the number of carbons. Activation energy is in the
* range (-1000 - -1500 K-1) and was set to the mean value -1250 (K-1). 
* Final expression is therefore : 
* k=3.43E-13*(1-exp(-0.245*nca))*exp(1250/T)
        arrh1 = 3.43E-13* ( 1-exp(-0.245*REAL(nca)) ) 
        arrh2 = 0. 
        arrh3 = -1250
      endif

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


* Generally, 1 channel (alkyl_peroxy to hydroperoxide)
* For case with a beta-carbonyl [RCOCR'R"(OO.)] add 2nd channel
* (alkyl_peroxy to alkoxy + OH, BR = 0.2) 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.

* first 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(1),nring)
      CALL stdchm(pchem(1))

      s(1) = 1.0
      s(2) = 0.0

* check if product is substituted alkene
* (and decompose it if necessary)
      !IF(INDEX(pchem(1),'CdH(OH)').NE.0) THEN
      IF(INDEX(pchem(1),'CdH(').NE.0.OR.
     &   INDEX(pchem(1),'Cd(').NE.0) THEN
        DO i=1,mca
          IF(INDEX(rdcktprod(1,i),' ').eq.1) THEN
            CALL alkcheck(pchem(1),rdcktprod(1,i))
            tempkcp=rdcktprod(1,i)
            IF(INDEX(pchem(1),'.').NE.0) THEN
              tempkc = pchem(1)
              CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
              IF (nip.EQ.2) THEN
                PRINT*,"error in ro2 + ho2, decompose subst. alkene"
                PRINT*,"------------> too many products!"
                STOP
              ENDIF
              pchem(1) = rdckprod(1)
            ENDIF
            CALL stdchm(pchem(1))
            s(1) = 1.
            rdcktprod(1,i) = tempkcp
            EXIT
          ENDIF
        ENDDO
* second product is coproduct from alkcheck
* (code will stop if radchk produces additional coproducts)
        np = np + 1
        s(np) = 1.
        p(np) = tempkcp
      ENDIF

! secondary channel
! following Wennberg et al., 2018
      IF (ro2ho2_fg .eq. 2) THEN
! s(2) = 0.1*n_betaOH+0.5*(n_alphaC=O + n_betaONO2 + n_betaOOH) 
!       +0.05*n_gammaOH + 0.25 *(n_betaC=O + n_gammaONO2 + n_gammaOOH)
! primary RO2 s(2) <= 0.2
! secondary RO2 s(2) <= 0.75
! tertiary RO2 s(2) <= 1.

! count functions in alpha and beta
! calculate chemmap for info about functinos
        CALL chemmap(chem,ngrp,group,bond,nfunc,nodetype,
     &            alifun,cdfun,arofun,mapfun,funflg,
     &            tabester,nfcd,nfcr,ierr)
        IF (ierr.ne. 0) THEN
          write(6,*) "error in ro2.f after chemmap call"
          write(6,*) "ierr=", ierr
          stop
        endif
        
! and use node positions with respect to peroxy radical
! calculated with abcde_map
        CALL abcde_map(bond, ip, ngrp, nabcde, tabcde)
        
        n_alphaCO = 0
        n_betaOH = 0
        n_betaONO2 = 0
        n_betaOOH = 0
        n_gammaOH = 0
        n_betaCO = 0
        n_gammaONO2 = 0
        n_gammaOOH = 0
        do j=1, nabcde(2)
          k = tabcde(2, j, 2)
          n_alphaCO = n_alphaCO + mapfun(k , 1, 9)  +
     &                mapfun(k, 1, 10)
          n_betaOH = n_betaOH + mapfun(k, 1, 1)
          n_betaONO2 = n_betaONO2 + mapfun(k, 1, 3)
          n_betaOOH = n_betaOOH + mapfun(k, 1, 4)
        enddo
        do j=1, nabcde(3)
          k = tabcde(3, j, 3)
          n_betaCO = n_betaCO + mapfun(k, 1, 9)  +
     &                mapfun(k, 1, 10)
          n_gammaOH = n_gammaOH + mapfun(k, 1, 1)
          n_gammaONO2 = n_gammaONO2 + mapfun(k, 1, 3)
          n_gammaOOH = n_gammaOOH + mapfun(k, 1, 4)
        enddo
                
! calculate s(2)
        s(2) = 0.1*n_betaOH + 0.5*(n_alphaCO + n_betaONO2 + n_betaOOH)
     & + 0.05*n_gammaOH + 0.25*(n_betaCO + n_gammaONO2 + n_gammaOOH)
     
! determine if ro2 is primary, secondary or tertiary
        if (itype .ge. 1 .and. itype .le. 3) then
          s(2) = min(0.2, s(2))
        else if (itype .ge. 4 .and. itype .le. 6) then
          s(2) = min(0.75, s(2))
        else if (itype .ge. 7 .and. itype .le. 8) then
          s(2) = min(1.00, s(2))
        endif
        s(1) = 1. - s(2)
        
* secondary channel:  change (OO.) to (O.); coproduct = OH
        pold = alkyl_peroxy
        pnew = alkoxy
        CALL swap(group(ip),pold,tgroup(ip),pnew)

*rebuild and rename:
        CALL rebond(tbond,tgroup,pchem(2),nring)
        CALL stdchm(pchem(2))        
        p(3) = 'HO    '
        s(3) = s(2)        
        
! or old method for secondary channel
      ELSEIF (iyl.eq.1) THEN
* secondary channel:  change (OO.) to (O.); coproduct = OH
        pold = alkyl_peroxy
        pnew = alkoxy
        CALL swap(group(ip),pold,tgroup(ip),pnew)

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

        s(1) = 0.8
        s(2) = 0.2
        p(3) = 'HO    '
        s(3) = s(2)
      ENDIF
      
* third product : remove the counters
*      s(3) = -1.
*      IF (itype.EQ.1) p(2) = 'XP1O2'
*      IF (itype.EQ.2) p(2) = 'XP2O2'
*      IF (itype.EQ.3) p(2) = 'XP3O2'
*      IF (itype.EQ.4) p(2) = 'XS1O2'
*      IF (itype.EQ.5) p(2) = 'XS2O2'
*      IF (itype.EQ.6) p(2) = 'XS3O2'
*      IF (itype.EQ.7) p(2) = 'XT1O2'
*      IF (itype.EQ.8) p(2) = 'XT2O2'

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

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

* to estimate the maximum production of ROOH (dbrch), we assumed :
* [NO]min = 10ppt = 2.46E8 molec.cm-3
* [HO2]max = 1E9 molec.cm-3
* maximum is at 298K 1atm
c      brtio = brch * f298*1E9 / 
c     & (f298*1E9 + 2.7E-12*exp(360./298.)*2.46E8)
c      brtio = brch 
      brtio = brch*s(1)!*0.5
      CALL bratio(pchem(1),brtio,p(1),
     &                  dbrch,dict,namlst,
     &                  nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &                  nfn,namfn,chemfn)
      brtio = brch*s(2)!*0.5

      IF (s(2).NE.0) THEN
      CALL bratio(pchem(2),brtio,p(2),
     &                  dbrch,dict,namlst,
     &                  nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &                  nfn,namfn,chemfn)
      ENDIF

* rien n'est prevu pour le coproduit PARC5 s'il est produit ici
      IF (rdtcopchem.GT.0.) THEN
        WRITE(6,*)'WARNING from ro2.f'
        WRITE(6,*) 'Hydroperoxyde replaced by lumped species'
        WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF      
      
      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 ro2 + ho2'


* -------------------
* reaction with NO3
* -------------------
      IF (wtflag.NE.0) print*,'rxn with NO3'

* rate constant :
* Recommended by Atkinson 99 for C2H5O2 (Biggs et al., J. Chem. Soc.
* Faraday Trans., 817, 1995 ; Ray et al., J. Phys. Chem., 5737, 1996).
* C2H5O2 rate constant is used for every RO2 
      arrh1 =  2.3E-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,' '))

* ihange (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(1) = rdckprod(1)
      CALL stdchm(pchem(1))
      s(1) = 1.
      IF (nip.EQ.2) THEN
        pchem(2) = rdckprod(2)
        CALL stdchm(pchem(2))
        s(2) = sc(2)
        s(1) = sc(1)
      ENDIF

* check if product is substituted alkene
* (and decompose it if necessary)
      IF(INDEX(pchem(1),'CdH(').NE.0.OR.
     &   INDEX(pchem(1),'Cd(').NE.0) THEN
        DO i=1,mca
          IF(INDEX(rdcktprod(1,i),' ').eq.1) THEN
            CALL alkcheck(pchem(1),rdcktprod(1,i))
            tempkcp=rdcktprod(1,i)
            IF(INDEX(pchem(1),'.').NE.0) THEN
              tempkc = pchem(1)
              CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
              IF (nip.EQ.2) THEN
                PRINT*,"error in ro2 + no3, decompose subst. alkene"
                PRINT*,"------------> too many products!"
                STOP
              ENDIF
              pchem(1) = rdckprod(1)
            ENDIF
            CALL stdchm(pchem(1))
            s(1) = 1.
            rdcktprod(1,i) = tempkcp
            EXIT
          ENDIF
        ENDDO
* second product is coproduct from alkcheck
* (code will stop if radchk produces additional coproducts)
        np = np + 1
        s(np) = 1.
        p(np) = tempkcp
      ENDIF

      coprod(:) = rdcktprod(1,:)
      coprod2(:) = rdcktprod(2,:)

      brtio = brch*0.01 
      CALL bratio(pchem(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      np = 1
      IF (nip.EQ.2) THEN
      np = np + 1
      CALL bratio(pchem(np),brtio,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.
      p(np) = 'NO2  '

* third product : remove the counters
*      s(3) = -1.
*      IF (itype.EQ.1) p(3) = 'XP1O2'
*      IF (itype.EQ.2) p(3) = 'XP2O2'
*      IF (itype.EQ.3) p(3) = 'XP3O2'
*      IF (itype.EQ.4) p(3) = 'XS1O2'
*      IF (itype.EQ.5) p(3) = 'XS2O2'
*      IF (itype.EQ.6) p(3) = 'XS3O2'
*      IF (itype.EQ.7) p(3) = 'XT1O2'
*      IF (itype.EQ.8) p(3) = 'XT2O2'

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

* reactant
      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 ro2 + no3 abstraction'


* -------------------------------
* RO2+RO2 reactions
* -------------------------------
      IF (wtflag.NE.0) print*,'rxn with RO2'

* species for which the "self" reaction rates 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. 
c      wf=0.001 ! weighting factor applied  to bratio for the RO2+RO2 recations
      wf=0.02 ! weighting factor applied  to bratio for the RO2+RO2 recations
c      IF (itype.eq.4.OR.itype.eq.7.OR.itype.eq.8) GOTO 900
      IF (itype.eq.4) GOTO 900

* find the various products linked to the species (alkoxy,
* alcohol, carbonyl)

      p_rad(1)  = ' '
      p_rad(2)  = ' '
      p_ol   = ' '
      p_onyl = ' '

      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(1) = rdckprod(1)
      s(1) = 1.
      IF (nip.EQ.2) THEN
        p_rad(2) = rdckprod(2)
        CALL stdchm(p_rad(2))
      ENDIF
      coprod(:) = rdcktprod(1,:)
      coprod2(:) = rdcktprod(2,:)
      CALL stdchm(p_rad(1))

      pold = alkyl_peroxy
      pnew = hydroxy
      CALL swap(group(ip),pold,tgroup(ip),pnew)
      CALL rebond(tbond,tgroup,p_ol,nring)
      CALL stdchm(p_ol)
! ric august 2015 : this structure appears in isoprene chemistry from MCM 
!      IF((INDEX(p_ol,'CdH(OH)').NE.0).OR.
!     &   (INDEX(p_ol,'Cd(OH)').NE.0)) THEN
!        CALL alkcheck(p_ol,prod_ol)
!      ENDIF
      
* check if product is substituted alkene
* (and decompose it if necessary)
*  >Cd=Cd(OH)-    -> >CH-CO-   (no call to radchk, no coproduct)
        IF(INDEX(p_ol,'CdH(').NE.0.OR.
     &     INDEX(p_ol,'Cd(').NE.0) THEN
          DO i=1,mca
            IF(INDEX(rdcktprod(1,i),' ').eq.1) THEN
              CALL alkcheck(p_ol,rdcktprod(1,i))
              p_ol = rdckprod(1)
              CALL stdchm(p_ol)
              EXIT
            ENDIF
          ENDDO
        ENDIF

      IF (INDEX(group(ip),'CH2').NE.0) THEN
         tgroup(ip)=aldehyde
         CALL rebond(tbond,tgroup,p_onyl,nring)
         CALL stdchm(p_onyl)
      ELSE IF (INDEX(group(ip),'CH').NE.0) THEN
         pold = 'CH'
         pnew= 'CO'
         CALL swap(group(ip),pold,tgroup(ip),pnew)
         pold = alkyl_peroxy
         pnew= ' '
         tempkg=tgroup(ip)
         CALL swap(tempkg,pold,tgroup(ip),pnew)
         CALL rebond(tbond,tgroup,p_onyl,nring)
         CALL stdchm(p_onyl)
      ENDIF

* RO2 = CH3(OO.)  (NAME = CH3O2)
* ===============================

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

* If radical channel for the processed RO2 is not equal to 1, then divide
* the molecular channel in two parts. However, if carbonyl formation
* is not allowed (e.g. in case of substitution at the peroxy position)
* then assume that reaction goes true radical channel. In what follows,
* rmol1 is for "H received" (from CH3O2) and rmol2 is for "H given" (to
* CH3O2)
      IF (ro2dat(itype,3).eq.1.) THEN
         rmol1=rmol
         rmol2=0.
      ELSE
         rmol1=rmol/2.       
         rmol2=rmol/2.
         IF (p_onyl(1:2).eq.'  ') THEN
           rrad=rrad+rmol2
           rmol2=0.
         ENDIF  
      ENDIF     

* change stoichiometric coefficient in such a way that sum = 1
* when writing the results. This is done by changing rrad.
      WRITE(wrtnum,'(f5.2)') rmol1 
      READ(wrtnum,'(f5.2)') rmol1
      WRITE(wrtnum,'(f5.2)') rmol2 
      READ(wrtnum,'(f5.2)') rmol2
      check=1.-rmol1-rmol2
      IF (abs((check-rrad)/rrad).gt.0.03) THEN
         WRITE(6,*) '--error-- in RO2 (treating CH3O2 reaction)'
         WRITE(6,*) 'something wrong in overwriting stoic. coef.'
         WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF
      rrad=check

* rate constant : (2*geometric mean)
      arrh1 =  2.*( (ch3o2dat(1)*ro2dat(itype,1)) )**0.5
      arrh2 =  0.
      arrh3 =  (ch3o2dat(2)+ro2dat(itype,2))/2. 
      IF (wtflag.NE.0) THEN
        write(6,*) 'arrh1=',arrh1
        write(6,*) 'arrh2=',arrh2
        write(6,*) 'arrh3=',arrh3
      ENDIF        

      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)
      r(2) = 'MEPERO'

* radical channel
      s(1)=rrad
c      brtio = brch * s(1) *f298*1E9/(f298*1E9 + 2.7E-12*exp(360./298.)*
c     &     2.46E8) 
      brtio=brch*wf
      CALL bratio(p_rad(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      np = 1
      IF (nip.EQ.2) THEN
        np = np +1
      CALL bratio(p_rad(2),brtio,p(2),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
        s(1)=rrad*sc(1)
        s(2)=rrad*sc(2)
      ENDIF
c      s(2)=rrad
c      p(2) = 'CH3O '

* "H received" from CH3O2 channel
      np = np + 1
      s(np)=rmol1
c      brtio = brch * s(2) *f298*1E9/(f298*1E9 + 2.7E-12*exp(360./298.)*
c     &     2.46E8) 
      IF ((itype.gt.3).and.(itype.lt.9))THEN
        brtio = brch*0.2*wf
      ELSE
        brtio = brch*0.5*wf
      ENDIF
      CALL bratio(p_ol,brtio,p(np),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

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

* "H given" to CH3O2 channel
c      np=2
      IF (rmol2.gt.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rmol2
c        brtio = brch * s(np)*f298*1E9/(f298*1E9 + 2.7E-12*exp(360./298.)
c     &     *2.46E8) 
        IF ((itype.gt.3).and.(itype.lt.9))THEN
          brtio = brch*0.2*wf
        ELSE
          brtio = brch*0.5*wf
        ENDIF
        CALL bratio(p_onyl,brtio,p(np),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)

c        CALL addprod(np,progname,rdct(lco+1:lcf))
c        s(np) = rmol2
c        p(np) = 'CH3OH'
      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
      IF (nip.EQ.2) THEN
        DO i=1,mca
          IF (coprod2(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(2)
            p(np) = coprod2(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
          ENDIF
        ENDDO
      ENDIF
      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
      IF (rdtcopchem3.GT.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rdtcopchem3
        p(np) = copchem3
      ENDIF

* last product : remove the counters
      CALL addprod(np,progname,rdct(lco+1:lcf))
*      s(np) = -1.
*      IF (itype.EQ.1) p(np) = 'XP1O2'
*      IF (itype.EQ.2) p(np) = 'XP2O2'
*      IF (itype.EQ.3) p(np) = 'XP3O2'
*      IF (itype.EQ.4) p(np) = 'XS1O2'
*      IF (itype.EQ.5) p(np) = 'XS2O2'
*      IF (itype.EQ.6) p(np) = 'XS3O2'
*      IF (itype.EQ.7) p(np) = 'XT1O2'
*      IF (itype.EQ.8) p(np) = 'XT2O2'

* 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
        IF (j.eq.4) GOTO 200
        IF (j.eq.7) GOTO 200
        IF (j.eq.8) GOTO 200

* compute the branching ratio
        rmol=( (1.-ro2dat(j,3)) + (1.-ro2dat(itype,3)) )/2.
        rrad=1.-rmol

* If radical channel for the processed RO2 is not equal to 1, then divide
* the molecular channel in two parts. However, if carbonyl formation
* is not allowed (e.g. in case of substitution at the peroxy position)
* then assume that reaction goes through radical channel. In what follows,
* rmol1 is for "H received" (from counter) and rmol2 is for "H given" (to
* counter)
        IF (rrad.gt.0.999) THEN
           rmol1=0.
           rmol2=0.
           rrad=1.
        ELSE IF (ro2dat(itype,3).eq.1.) THEN
           rmol1=rmol
           rmol2=0.
        ELSE IF (ro2dat(j,3).eq.1.) THEN
           rmol1=0.
           rmol2=1.
        ELSE
           rmol1=rmol/2.       
           rmol2=rmol/2.
           IF (p_onyl(1:2).eq.'  ') THEN
             rrad=rrad+rmol2
             rmol2=0.
           ENDIF  
        ENDIF     

* change stoichimetric coefficients in such a way that sum = 1
* when writing the results. This is done by changing rrad.
      WRITE(wrtnum,'(f5.2)') rmol1 
      READ(wrtnum,'(f5.2)') rmol1
      WRITE(wrtnum,'(f5.2)') rmol2 
      READ(wrtnum,'(f5.2)') rmol2
      check=1.-rmol1-rmol2
      IF (abs((check-rrad)/rrad).gt.0.03) THEN
         WRITE(6,*) '--error-- in RO2 (treating RO2 reaction)'
         WRITE(6,*) 'something wrong in overwriting stoic. coef.'
         WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF
      rrad=check

* rate constant : (2*geometric mean)
        arrh1 =  2.*( (ro2dat(j,1)*ro2dat(itype,1)) )**0.5
        arrh2 =  0.
        arrh3 =  ( ro2dat(j,2)+ro2dat(itype,2) )/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)='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)='PERO'
        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/(f298*1E9 + 2.7E-12*exp(360./298.)*
c     &     2.46E8) 
        brtio = brch*wf 
        CALL bratio(p_rad(1),brtio,p(1),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
        np = 1
        IF (nip.EQ.2) THEN
          np = np +1
        CALL bratio(p_rad(2),brtio,p(2),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
          s(1)=rrad*sc(1)
          s(2)=rrad*sc(2)
        ENDIF

* "H received from counter" channel
c        np=1
        IF (rmol1.gt.0.) THEN
          CALL addprod(np,progname,rdct(lco+1:lcf))
          s(np) = rmol1
c          brtio = brch * s(np)*f298*1E9/(f298*1E9 + 2.7E-12
c     &    *exp(360./298.)* 2.46E8) 
          IF ((itype.gt.3).and.(itype.lt.9))THEN
            brtio = brch*0.2*wf
          ELSE
            brtio = brch*0.5*wf
          ENDIF
          CALL bratio(p_ol,brtio,p(np),
     &              dbrch,dict,namlst,
     &              nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &              nfn,namfn,chemfn)
        ENDIF
        IF (rdtcopchem.GT.0.) THEN
          CALL addprod(np,progname,rdct(lco+1:lcf))
          s(np) = rdtcopchem
          p(np) = copchem
        ENDIF

* "H given to counter" channel
        IF (rmol2.gt.0.) THEN
          CALL addprod(np,progname,rdct(lco+1:lcf))
          s(np) = rmol2
c          brtio = brch * s(np)*f298*1E9/(f298*1E9 + 2.7E-12*
c     &    exp(360./298.)* 2.46E8) 
          IF ((itype.gt.3).and.(itype.lt.9))THEN
            brtio = brch*0.2*wf
          ELSE
            brtio = brch*0.5*wf
          ENDIF
          CALL bratio(p_onyl,brtio,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
        IF (nip.EQ.2) THEN
        DO i=1,mca
           IF (coprod2(i)(1:1).NE.' ') THEN
              CALL addprod(np,progname,rdct(lco+1:lcf))
              s(np) = s(2)
              p(np) = coprod2(i)
              IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
           ENDIF
        ENDDO
        ENDIF
        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)
        idreac=0
c        idreac=2
        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 ro2 reaction with counter 
     &        number',j

200   CONTINUE

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

* compute the branching ratio
      IF (ro2dat(itype,3).eq.1.) THEN
         rrad=1.
         rmol1=0.
         rmol2=0.
      ELSE
         rmol1=0.       
         rmol2=( 1.-ro2dat(itype,3) )/2.
         rrad=1.-rmol2
         IF (p_onyl(1:2).eq.'  ') THEN
           rrad=1.
           rmol2=0.
         ENDIF  
      ENDIF     

* change stoichimetric coefficient in such a way that sum make 1
* when writing the results. This is done by changing rrad.
      WRITE(wrtnum,'(f5.2)') rmol1 
      READ(wrtnum,'(f5.2)') rmol1
      WRITE(wrtnum,'(f5.2)') rmol2 
      READ(wrtnum,'(f5.2)') rmol2
      check=1.-rmol1-rmol2
      IF (abs((check-rrad)/rrad).gt.0.03) THEN
         WRITE(6,*) '--error-- in RO2 (treating RCO3 reaction)'
         WRITE(6,*) 'something wrong in overwriting stoic. coef.'
         WRITE(99,*) 'ro2',rdct(lco+1:lcf) !STOP
      ENDIF
      rrad=check


* rate constant : Data available shows that the reaction rate
* is close to 1E-11, independent of the RO2 being considered. This
* rate constant is used, instead of the classical geometrical mean
C      arrh1 =  2.*( (rco3dat(1)*ro2dat(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)
      r(2) = 'PERO9'
c      nlabel=409

* radical channel
      s(1)=rrad
c      brtio = brch * s(1) *f298*1E9/(f298*1E9 + 2.7E-12
c     &    *exp(360./298.)* 2.46E8) 
      brtio = brch*wf
      CALL bratio(p_rad(1),brtio,p(1),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
      np = 1
      IF (nip.EQ.2) THEN
        np = np +1
      CALL bratio(p_rad(2),brtio,p(2),
     &            dbrch,dict,namlst,
     &            nhldvoc,holdvoc,nhldrad,holdrad,level,stabl,
     &            nfn,namfn,chemfn)
        s(1)=rrad*sc(1)
        s(2)=rrad*sc(2)
      ENDIF

* "H given" to RCO3
c      np=1
      IF (rmol2.gt.0.) THEN
        CALL addprod(np,progname,rdct(lco+1:lcf))
        s(np) = rmol2
c        brtio = brch * s(np) *f298*1E9/(f298*1E9 + 2.7E-12
c     &    *exp(360./298.)* 2.46E8) 
        IF ((itype.gt.3).and.(itype.lt.9))THEN
          brtio = brch*0.2*wf
        ELSE
          brtio = brch*0.5*wf
        ENDIF
        CALL bratio(p_onyl,brtio,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
      IF (nip.EQ.2) THEN
      DO i=1,mca
         IF (coprod2(i)(1:1).NE.' ') THEN
            CALL addprod(np,progname,rdct(lco+1:lcf))
            s(np) = s(2)
            p(np) = coprod2(i)
            IF (wtopeflag.EQ.1) write(10,'(23X,A1,A6)')'G',p(np)
         ENDIF
      ENDDO
      ENDIF
      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
900   CONTINUE
99    CONTINUE
* end RO2
      RETURN
      END
