Index: trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2164)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2170)
@@ -26,5 +26,5 @@
       use conc_mod, only: mmean ! mean molecular mass of the atmosphere
       use comcstfi_h, only: pi
-      use photolysis_mod, only: jonline, init_photolysis
+      use photolysis_mod, only: init_photolysis, nphot
       use iono_h, only: temp_elect
 
@@ -38,14 +38,4 @@
 !  Prepare the call for the photochemical module, and send back the
 !  tendencies from photochemistry in the chemical species mass mixing ratios
-!
-!   Author:   Sebastien Lebonnois (08/11/2002)
-!   -------
-!    update 12/06/2003 for water ice clouds and compatibility with dust
-!    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
-!    update 03/05/2005 cosmetic changes (Franck Lefevre)
-!    update sept. 2008 identify tracers by their names (Ehouarn Millour)
-!    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
-!    update 16/03/2012 optimization (Franck Lefevre)
-!    update 11/12/2013 optimization (Franck Lefevre)
 !
 !   Arguments:
@@ -155,4 +145,11 @@
       integer :: ig_vl1
 
+      integer :: nb_reaction_3_max     ! number of quadratic reactions
+      integer :: nb_reaction_4_max     ! number of bimolecular reactions
+      integer :: nquench               ! number of quenching + heterogeneous reactions
+      integer :: nphotion              ! number of photoionizations
+      integer :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
+
+
       real    :: latvl1, lonvl1
       real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
@@ -164,6 +161,11 @@
       real    :: kb                    ! boltzmann constant
 
-      logical,save :: firstcall = .true.
-      logical,save :: depos = .false.  ! switch for dry deposition
+      logical, save :: firstcall = .true.
+      logical, save :: depos       ! switch for dry deposition
+      logical, save :: ionchem     ! switch for ionospheric chemistry
+      logical, save :: jonline     ! switch for online photodissociation rates or lookup table
+      logical, save :: unichim     ! only one unified chemical scheme at all 
+                                   ! layers (default), or upper atmospheric 
+                                   ! scheme in the thermosphere
 
 !     for each column of atmosphere:
@@ -181,14 +183,13 @@
       real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
       real :: jh2o(nlayer)         ! Photodissociation rate H2O->H+OH (s-1)
-      real :: em_no(nlayer)        !  NO nightglow emission rate
-      real :: em_o2(nlayer)        !  O2 nightglow emission rate      
-
-      integer :: iter(nlayer)      !  Number of chemical iterations
-                                   !  within one physical timestep 
+      real :: em_no(nlayer)        ! NO nightglow emission rate
+      real :: em_o2(nlayer)        ! O2 nightglow emission rate      
+
+      integer :: iter(nlayer)      ! Number of chemical iterations
+                                   ! within one physical timestep 
 
 !     for output:
 
       logical :: output             ! to issue calls to writediagfi and stats
-      parameter (output = .true.)
       real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
       real :: jh2o_3d(ngrid,nlayer)  ! Photodissociation rate H2O->H+OH (s-1)
@@ -198,16 +199,17 @@
                                     !  within one physical timestep
 
-      logical :: unichim            ! only one, unified chemical scheme at all 
-                                    ! layers (default), or upper atmospheric 
-                                    ! scheme in the thermosphere?
+!=======================================================================
+!     main dashboard for the chemistry
+!=======================================================================
+
+      unichim = .true.     ! true : unified chemistry ! false : separate models in lower and upper atmosphere
+      jonline = .true.     ! true : on-line calculation of photodissociation rates ! false : lookup table
+      ionchem = .false.    ! switch for ionospheric chemistry
+      depos   = .false.    ! switch for dry deposition
+      output  = .false.    ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)
 
 !=======================================================================
 !     initialization of the chemistry (first call only)
 !=======================================================================
-
-      !call only unified chemistry (default), or also upper atmospheric model
-      !unichim = .true.   unified chemistry
-      !unichim = .false.  2 different models
-      unichim = .true.
 
       if (firstcall) then
@@ -388,91 +390,175 @@
          i_o2plus = igcm_o2plus
          if (i_o2plus == 0) then
-            write(*,*) "calchim: Error, no O2+ tracer !!!"
-            stop
+            write(*,*) "calchim: no O2+ tracer "
+            write(*,*) "Only neutral chemistry"
          else
             nbq = nbq + 1
             niq(nbq) = i_o2plus
+            ionchem = .true.
+            write(*,*) "calchim: O2+ tracer found in traceur.def"
+            write(*,*) "Ion chemistry included"
          end if
          i_co2plus = igcm_co2plus
-         if (i_co2plus == 0) then
-            write(*,*) "calchim: Error, no CO2+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_co2plus
-         end if
+         if(ionchem) then
+            if (i_co2plus == 0) then
+               write(*,*) "calchim: Error, no CO2+ tracer !!!"
+               write(*,*) "CO2+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_co2plus
+            end if
+         else
+            if (i_co2plus /= 0) then
+               write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+               stop
+            endif
+         endif
          i_oplus=igcm_oplus
-         if (i_oplus == 0) then
-            write(*,*) "calchim: Error, no O+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_oplus
-         end if
+         if(ionchem) then
+            if (i_oplus == 0) then
+               write(*,*) "calchim: Error, no O+ tracer !!!"
+               write(*,*) "O+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_oplus
+            end if
+         else
+            if (i_oplus /= 0) then
+               write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+               stop
+            endif
+         endif
          i_noplus=igcm_noplus
-         if (i_noplus == 0) then
-            write(*,*) "calchim: Error, no NO+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_noplus
-         end if
+         if(ionchem) then
+            if (i_noplus == 0) then
+               write(*,*) "calchim: Error, no NO+ tracer !!!"
+               write(*,*) "NO+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_noplus
+            end if
+         else
+            if (i_noplus /= 0) then
+               write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_coplus=igcm_coplus
-         if (i_coplus == 0) then
-            write(*,*) "calchim: Error, no CO+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_coplus
-         end if
+         if(ionchem) then
+            if (i_coplus == 0) then
+               write(*,*) "calchim: Error, no CO+ tracer !!!"
+               write(*,*) "CO+ is needed if O2+ is in traceur.def"
+               stop               
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_coplus
+            end if
+         else
+            if (i_coplus /= 0) then
+               write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_cplus=igcm_cplus
-         if (i_cplus == 0) then
-            write(*,*) "calchim: Error, no C+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_cplus
-         end if
+         if(ionchem) then
+            if (i_cplus == 0) then
+               write(*,*) "calchim: Error, no C+ tracer !!!"
+               write(*,*) "C+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_cplus
+            end if
+         else
+            if (i_cplus /= 0) then
+               write(*,*) "calchim: Error: C+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_n2plus=igcm_n2plus
-         if (i_n2plus == 0) then
-            write(*,*) "calchim: Error, no N2+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n2plus
-         end if
+         if(ionchem) then
+            if (i_n2plus == 0) then
+               write(*,*) "calchim: Error, no N2+ tracer !!!"
+               write(*,*) "N2+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_n2plus
+            end if
+         else
+            if (i_n2plus /= 0) then
+               write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_nplus=igcm_nplus
-         if (i_nplus == 0) then
-            write(*,*) "calchim: Error, no N+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_nplus
-         end if
+         if(ionchem) then
+            if (i_nplus == 0) then
+               write(*,*) "calchim: Error, no N+ tracer !!!"
+               write(*,*) "N+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_nplus
+            end if
+         else
+            if (i_nplus /= 0) then
+               write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_hplus=igcm_hplus
-         if (i_hplus == 0) then
-            write(*,*) "calchim: Error, no H+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_hplus
-         end if
+         if(ionchem) then
+            if (i_hplus == 0) then
+               write(*,*) "calchim: Error, no H+ tracer !!!"
+               write(*,*) "H+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_hplus
+            end if
+         else
+            if (i_hplus /= 0) then
+               write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_hco2plus=igcm_hco2plus
-         if (i_hco2plus == 0) then
-            write(*,*) "calchim: Error, no HCO2+ tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_hco2plus
-         end if
+         if(ionchem) then
+            if (i_hco2plus == 0) then
+               write(*,*) "calchim: Error, no HCO2+ tracer !!!"
+               write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_hco2plus
+            end if
+         else
+            if (i_hco2plus /= 0) then
+               write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          i_elec = igcm_elec
-         if (i_elec == 0) then
-            write(*,*) "calchim: Error, no e- tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_elec
-         end if
-         
+         if(ionchem) then
+            if (i_elec == 0) then
+               write(*,*) "calchim: Error, no e- tracer !!!"
+               write(*,*) "e- is needed if O2+ is in traceur.def"
+               stop
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_elec
+            end if
+         else
+            if (i_elec /= 0) then
+               write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
          write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
                
@@ -523,6 +609,7 @@
 
 !           search for switch index between regions  
-            if(unichim) then
-               lswitch=nlayer+1
+
+            if (unichim) then
+               lswitch = nlayer + 1
             else
                if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
@@ -540,11 +627,30 @@
 !=======================================================================
 
-!        chemistry: only one scheme at all layers
-         
          if (photochem) then
-            call photochemistry(nlayer,nq,                            &
-                           ig,lswitch,zycol,szacol,ptimestep,         &
-                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,   &
-                           dist_sol,zday,surfdust1d,surfice1d,        &
+            ! set number of reactions, depending on ion chemistry or not
+            if (ionchem) then
+               nb_reaction_4_max = 60   ! set number of bimolecular reactions
+               nb_reaction_3_max = 6    ! set number of quadratic reactions
+               nquench           = 9    ! set number of quenching + heterogeneous reactions
+               nphotion          = 18   ! set number of photoionizations
+            else
+               nb_reaction_4_max = 31   ! set number of bimolecular reactions
+               nb_reaction_3_max = 6    ! set number of quadratic reactions
+               nquench           = 9    ! set number of quenching + heterogeneous reactions
+               nphotion          = 0    ! set number of photoionizations
+            end if
+
+!        nb_phot_max is the total number of processes that are treated
+!        numerically as a photolysis:
+
+            nb_phot_max = nphot + nphotion + nquench
+
+!        call main photochemistry routine
+
+            call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max,  &
+                           nb_reaction_4_max, nb_phot_max, nphotion,      &
+                           jonline, ig,lswitch,zycol,szacol,ptimestep,    &
+                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
+                           dist_sol,zday,surfdust1d,surfice1d,            &
                            jo3,jh2o,taucol,iter)
 
@@ -556,4 +662,5 @@
                iter_3d(ig,l) = iter(l)
             end do
+
 !        condensation of h2o2
 
@@ -562,15 +669,18 @@
                          ztemp,zycol,dqcloud,dqscloud)
 
-            if(.not.unichim) then
-               chemthermod=3   !C/O/H/N/ions
+!        case of separate photochemical model in the thermosphere
+
+            if (.not.unichim) then
+               chemthermod = 3   !C/O/H/N/ions
                call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
                     zdens,zpress,zlocal,szacol,ptimestep,zday,&
                     em_no,em_o2)
-               do l=1,nlayer
-                  emission_no(ig,l)=em_no(l)
-                  emission_o2(ig,l)=em_o2(l)
-               enddo
-            end if
-         end if
+               do l = 1,nlayer
+                  emission_no(ig,l) = em_no(l)
+                  emission_o2(ig,l) = em_o2(l)
+               end do
+            end if
+
+         end if  ! photochem
 
 !        dry deposition
@@ -627,5 +737,6 @@
             call writediagfi(ngrid,'iter','iterations',  &
                              ' ',3,iter_3d(1,1))
-            if(.not.unichim) then
+
+            if (.not. unichim) then
                call writediagfi(ngrid,'emission_no',        &
                     'NO nightglow emission rate','cm-3 s-1',3,emission_no)
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2164)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2170)
@@ -13,14 +13,10 @@
 !*****************************************************************
 
-subroutine photochemistry(nlayer, nq,                                   &
-                          ig, lswitch, zycol, sza, ptimestep, press,    &
-                          alt, temp, temp_elect,dens, zmmean,           &
-                          dist_sol, zday,                               &
+subroutine photochemistry(nlayer, nq, nesp, ionchem, nb_reaction_3_max,        &
+                          nb_reaction_4_max, nb_phot_max, nphotion,            &
+                          jonline, ig, lswitch, zycol, sza, ptimestep, press,  &
+                          alt, temp, temp_elect, dens, zmmean,                 &
+                          dist_sol, zday,                                      &
                           surfdust1d, surfice1d, jo3, jh2o,tau, iter)
-
-use photolysis_mod, only : nb_phot_max,       &
-                           nb_reaction_3_max, &
-                           nb_reaction_4_max, &
-                           jonline
 
 use param_v4_h, only: jion
@@ -35,5 +31,16 @@
 
 integer, intent(in) :: nlayer ! number of atmospheric layers
-integer, intent(in) :: nq     ! number of tracers
+integer, intent(in) :: nq     ! number of tracers in traceur.def
+integer, intent(in) :: nesp   ! number of traceurs in chemistry
+integer, intent(in) :: nb_reaction_3_max   
+                              ! number of quadratic reactions
+integer, intent(in) :: nb_reaction_4_max
+                              ! number of bimolecular reactions
+integer, intent(in) :: nb_phot_max
+                              ! number of reactions treated numerically as photodissociations
+integer, intent(in) :: nphotion
+                              ! number of photoionizations
+logical, intent(in) :: ionchem! switch for ion chemistry
+logical, intent(in) :: jonline! switch for on-line calculation of photolysis rates
 integer :: ig                 ! grid point index
       
@@ -69,6 +76,4 @@
 !     local:
 !===================================================================
-
-integer, parameter :: nesp = 28  ! number of species in the chemical code
 
 integer :: phychemrat            ! (physical timestep)/(nominal chemical timestep)
@@ -109,5 +114,4 @@
 integer,parameter :: i_elec    = 28
 
-
 integer :: ilay
 
@@ -145,8 +149,9 @@
 if (firstcall) then
    print*,'photochemistry: initialize indexes'
-   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
-               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
+   call indice(nb_reaction_3_max,nb_reaction_4_max,                 &
+               nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
+               i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
+               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,            &
+               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
                i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
    firstcall = .false.
@@ -157,5 +162,5 @@
 !===================================================================
 
-call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
+call gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,      &
                i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
                i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
@@ -169,9 +174,9 @@
 !===================================================================
 
-jionos= .true.
+jionos = .true.
 
 if (jonline) then
    if (sza <= 113.) then ! day at 300 km
-      call photolysis_online(nlayer, alt, press, temp, zmmean,          &
+      call photolysis_online(nlayer, nb_phot_max, alt, press, temp, zmmean, &
                              i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
                              i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
@@ -179,5 +184,5 @@
                              tau, sza, dist_sol, v_phot)
 
-      if(jionos) then
+      if (jionos .and. ionchem) then
          call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
          do ilay=1,lswitch-1
@@ -246,6 +251,6 @@
 else
    tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
-   call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & 
-                   rm(:,i_co2), rm(:,i_o3), v_phot)
+   call photolysis(nlayer, nb_phot_max, lswitch, press, temp, sza, tau,  &
+                   zmmean, dist_sol,rm(:,i_co2), rm(:,i_o3), v_phot)
 end if
 ! save o3 and h2o photolysis for output
@@ -267,7 +272,8 @@
 hetero_ice  = .true.
 
-call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2),   &
-                   c(:,i_o), c(:,i_n2), press, temp, temp_elect,      &
-                   hetero_dust, hetero_ice,                           &
+call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
+                   nb_phot_max, nphotion, lswitch, dens, c(:,i_co2),      &
+                   c(:,i_o2), c(:,i_o), c(:,i_n2), press, temp,           &
+                   temp_elect, hetero_dust, hetero_ice,                   &
                    surfdust1d, surfice1d, v_phot, v_3, v_4)
 
@@ -307,5 +313,7 @@
 !  first-guess: fill matrix
 
-   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
+   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer,              &
+                    nb_reaction_3_max, nb_reaction_4_max, nb_phot_max,    &
+                    v_phot, v_3, v_4)
 
 !  adaptative evaluation of the sub time step
@@ -349,15 +357,19 @@
    cold(:)   = c(ilev,:)
    c(ilev,:) = cnew(:)
-   !Mod FGG, July 2019
-   !Force charge neutrality
-   if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
-        c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+ &
-        c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)) then
-      c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
-           c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+ &
-           c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)
-!      write(*,*)'photochemistry/359'
-!      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
-   endif
+
+!  force charge neutrality (mod fgg, july 2019)
+
+   if (ionchem) then
+      if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+&
+           c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
+           c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)) then
+         c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
+              c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+              &
+              c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+               &
+              c(ilev,i_hco2plus)
+         !      write(*,*)'photochemistry/359'
+         !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
+      end if
+   end if
    cnew(:)   = 0.
 
@@ -375,10 +387,10 @@
 !===================================================================
 
-call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
-               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,      &
-               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,   &
-               i_n2plus, i_nplus, i_hplus, i_hco2plus,        &
+call chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp,      &
+               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
+               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
+               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
+               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
+               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
                i_elec, dens, c) 
 contains
@@ -463,5 +475,5 @@
 
 #ifdef LAPACK
-      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
+   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
 #else
    write(*,*) "photochemistry error, missing LAPACK routine dgesv"
@@ -505,10 +517,8 @@
 !======================================================================
 
- subroutine reactionrates(nlayer,                                  &
-                          lswitch, dens, co2, o2, o, n2, press, &
-                          t, t_elect,                              &
-                          hetero_dust, hetero_ice,                 &
-                          surfdust1d, surfice1d,                   &
-                          v_phot, v_3, v_4)
+ subroutine reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
+                          nb_phot_max, nphotion, lswitch, dens, co2, o2, o,      &
+                          n2, press, t, t_elect, hetero_dust, hetero_ice,        &
+                          surfdust1d, surfice1d, v_phot, v_3, v_4)
  
 !================================================================
@@ -524,7 +534,5 @@
 
 use comcstfi_h
-use photolysis_mod, only : nphot, nphotion, nb_phot_max, &
-                           nb_reaction_3_max,  &
-                           nb_reaction_4_max
+use photolysis_mod, only : nphot
 
 implicit none
@@ -535,4 +543,9 @@
 
 integer, intent(in)     :: nlayer            ! number of atmospheric layers
+integer, intent(in)     :: nb_reaction_3_max ! number of quadratic reactions
+integer, intent(in)     :: nb_reaction_4_max ! number of bimolecular reactions
+integer, intent(in)     :: nb_phot_max       ! number of reactions treated numerically as photodissociations
+integer, intent(in)     :: nphotion          ! number of photoionizations
+logical, intent(in)     :: ionchem
 integer                 :: lswitch           ! interface level between lower
                                              ! atmosphere and thermosphere chemistries
@@ -563,5 +576,5 @@
 integer          :: ilev
 integer          :: nb_phot, nb_reaction_3, nb_reaction_4
-real :: ak0, ak1, xpo, rate
+real :: ak0, ak1, xpo, rate, rate1, rate2
 real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
 real, dimension(nlayer) :: deq
@@ -586,5 +599,5 @@
 !----------------------------------------------------------------------
 
-      nb_phot       = nphot + nphotion ! initialised to the number of photolysis and photoionization rates
+      nb_phot       = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
       nb_reaction_3 = 0
       nb_reaction_4 = 0
@@ -1022,17 +1035,21 @@
 !     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
 
-!     jpl 2006
-
-!     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
-!     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
-!     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
-!     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
-
-!     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
-!     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
-!     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
-!     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
-
-!     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
+!     jpl 2015
+
+!     do ilev = 1,lswitch-1
+
+!        branch 1 : oh + co -> h + co2
+
+!        rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
+
+!        branch 2 : oh + co + m -> hoco + m
+
+!        ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
+!        ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
+!        rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
+!        xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
+
+!        e001(ilev) = rate1 + rate2*0.6**xpo
+!     end do
 
 !     joshi et al., 2006
@@ -1053,5 +1070,4 @@
          k1a = k1a0*((1. + y)/(1. + x))*fx
          k1b = k1b0*(1./(1.+x))*fx
-
          e001(ilev) = k1a + k1b
       end do
@@ -1069,266 +1085,275 @@
       v_4(:,nb_reaction_4) = e002(:)
 
-!---  i001: co2+ + o2 -> o2+ + co2
-
-!     aninich, j. phys. chem. ref. data 1993
-
-      i001(:) = 5.5e-11*(300./t_elect(:))**0.5
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i001(:)
-
-!---  i002: co2+ + o -> o+ + co2
-
-!     UMIST database
-
-      i002(:) = 9.6e-11
+!----------------------------------------------------------------------
+!     ionospheric reactions
+!     only if ionchem=true
+!----------------------------------------------------------------------
+
+      if (ionchem) then
+
+!---     i001: co2+ + o2 -> o2+ + co2
+
+!        aninich, j. phys. chem. ref. data 1993
+
+         i001(:) = 5.5e-11*(300./t_elect(:))**0.5
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i001(:)
+
+!---     i002: co2+ + o -> o+ + co2
+
+!        UMIST database
+
+         i002(:) = 9.6e-11
       
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i002(:)
-
-!---  i003: co2+ + o -> o2+ + co
-
-!     UMIST database
-
-      i003(:) = 1.64e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i003(:)
-
-!---  i004: o2+ + e- -> o + o
-
-!     Alge et al., J. Phys. B, At. Mol. Phys. 1983
-
-      i004(:) = 2.0e-7*(300./t_elect(:))**0.7
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i004(:)
-
-!---  i005: o+ + co2 -> o2+ + co
-
-!     UMIST database
-
-      i005(:) = 9.4e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i005(:)
-
-
-!---  i006: co2+ + e- -> co + o
-
-!     UMIST database
-
-      i006(:) = 3.8e-7*(300./t_elect(:))**0.5
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i006(:)
-
-
-!---  i007: co2+ + no -> no+ + co2
-
-!     UMIST database
-
-      i007(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i007(:)
-
-!---  i008: o2+ + no -> no+ + o2
-
-!     UMIST database
-
-      i008(:) = 4.6e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i008(:)
-
-!---  i009: o2+ + n2 -> no+ + no
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i002(:)
+
+!---     i003: co2+ + o -> o2+ + co
+
+!        UMIST database
+
+         i003(:) = 1.64e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i003(:)
+
+!---     i004: o2+ + e- -> o + o
+
+!        Alge et al., J. Phys. B, At. Mol. Phys. 1983
+
+         i004(:) = 2.0e-7*(300./t_elect(:))**0.7
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i004(:)
+
+!---     i005: o+ + co2 -> o2+ + co
+
+!        UMIST database
+
+         i005(:) = 9.4e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i005(:)
+
+
+!---     i006: co2+ + e- -> co + o
+
+!        UMIST database
+
+         i006(:) = 3.8e-7*(300./t_elect(:))**0.5
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i006(:)
+
+
+!---     i007: co2+ + no -> no+ + co2
+
+!        UMIST database
+
+         i007(:) = 1.2e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i007(:)
+
+!---     i008: o2+ + no -> no+ + o2
+
+!        UMIST database
+
+         i008(:) = 4.6e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i008(:)
+
+!---     i009: o2+ + n2 -> no+ + no
       
-!     Fox & Sung 2001
-
-      i009(:) = 1.0e-15
+!        Fox & Sung 2001
+
+         i009(:) = 1.0e-15
       
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i009(:)
-
-!---  i010: o2+ + n -> no+ + o
-
-!     Fox & Sung 2001
-
-      i010(:) = 1.0e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i010(:)
-
-!---  i011: o+ + n2 -> no+ + n
-
-!     Fox & Sung 2001
-
-      i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i011(:)
-
-!---  i012: no+ + e -> n + o
-
-!     UMIST database
-
-      i012(:) = 4.3e-7*(300./t_elect(:))**0.37
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i012(:)
-
-
-!---  i013: co+ + co2 -> co2+ + co
-
-!     UMIST database
-
-      i013(:) = 1.0e-9
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i013(:)
-
-
-!---  i014: co+ + o -> o+ + co
-
-!     UMIST database
-
-      i014(:) = 1.4e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i014(:)
-
-!---  i015: c+ + co2 -> co+ + co
-
-!     UMIST database
-
-      i015(:) = 1.1e-9
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i015(:)
-
-
-!---  i016: N2+ + co2 -> co2+ + N2
-
-!     Fox & Song 2001
-
-      i016(:) = 9.0e-10*(300./t_elect(:))**0.23
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i016(:)
-
-
-!---  i017: N2+ + o -> no+ + N
-
-!     Fox & Song 2001
-
-      i017(:) = 1.33e-10*(300./t_elect(:))**0.44
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i017(:)
-
-!---  i018: N2+ + co -> co+ + N2
-
-!     UMIST
-
-      i018(:) = 7.4e-11
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i018(:)
-
-!---  i019: N2+ + e -> N + N
-
-!     UMIST
-
-      i019(:) = 7.7e-7*(300./t_elect(:))**0.3
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i016(:)
-
-!---  i020: N2+ + o -> o+ + N2
-
-!     Fox & Song 2001
-
-      i020(:) = 7.0e-12*(300./t_elect(:))**0.23
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i020(:)
-
-!---  i021: N+ + co2 -> co2+ + N
-
-!     UMIST
-
-      i021(:) = 7.5e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i021(:)
-
-!---  i022: CO+ + H -> H+ + CO
-
-!     Fox & Sung 2001
-
-      i022(:) = 4.0e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i022(:)
-
-!---  i023: O+ + H -> H+ + O
-
-!     UMIST
-
-      i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i023(:)
-
-!---  i024: H+ + O -> O+ + H
-
-!     UMIST
-
-      i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i024(:)
-
-!---  i025: CO+ + H2 -> HCO2+ + H
-
-!     UMIST
-
-      i025(:) = 9.5e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i025(:)
-
-!---  i026: HCO2+ + e -> H + CO2
-
-!     UMIST
-
-      i026(:) = 1.75e-8*((300./t_elect(:))**0.5)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i026(:)
-
-!---  i027+i028: HCO2+ + e -> H + O + CO
-
-!     UMIST
-      !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
-      !i028: 0.5 (HCO2+ + e-) -> O + CO
-
-      i027(:) = 2.38e-7*((300./t_elect(:))**0.5)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i027(:)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i027(:)
-
-!---  i029: HCO2+ + e -> H + CO2
-
-!     UMIST
-
-      i029(:) = 9.45e-8*((300./t_elect(:))**0.5)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = i029(:)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i009(:)
+
+!---     i010: o2+ + n -> no+ + o
+
+!        Fox & Sung 2001
+
+         i010(:) = 1.0e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i010(:)
+
+!---     i011: o+ + n2 -> no+ + n
+
+!        Fox & Sung 2001
+
+         i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i011(:)
+
+!---     i012: no+ + e -> n + o
+
+!        UMIST database
+
+         i012(:) = 4.3e-7*(300./t_elect(:))**0.37
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i012(:)
+
+
+!---     i013: co+ + co2 -> co2+ + co
+
+!        UMIST database
+
+         i013(:) = 1.0e-9
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i013(:)
+
+
+!---     i014: co+ + o -> o+ + co
+
+!        UMIST database
+
+         i014(:) = 1.4e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i014(:)
+
+!---     i015: c+ + co2 -> co+ + co
+
+!        UMIST database
+
+         i015(:) = 1.1e-9
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i015(:)
+
+
+!---     i016: N2+ + co2 -> co2+ + N2
+
+!        Fox & Song 2001
+
+         i016(:) = 9.0e-10*(300./t_elect(:))**0.23
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i016(:)
+
+
+!---     i017: N2+ + o -> no+ + N
+
+!        Fox & Song 2001
+
+         i017(:) = 1.33e-10*(300./t_elect(:))**0.44
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i017(:)
+
+!---     i018: N2+ + co -> co+ + N2
+
+!        UMIST
+
+         i018(:) = 7.4e-11
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i018(:)
+
+!---     i019: N2+ + e -> N + N
+
+!        UMIST
+
+         i019(:) = 7.7e-7*(300./t_elect(:))**0.3
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i016(:)
+
+!---     i020: N2+ + o -> o+ + N2
+
+!        Fox & Song 2001
+
+         i020(:) = 7.0e-12*(300./t_elect(:))**0.23
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i020(:)
+
+!---     i021: N+ + co2 -> co2+ + N
+
+!        UMIST
+
+         i021(:) = 7.5e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i021(:)
+
+!---     i022: CO+ + H -> H+ + CO
+
+!        Fox & Sung 2001
+
+         i022(:) = 4.0e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i022(:)
+
+!---     i023: O+ + H -> H+ + O
+
+!        UMIST
+
+         i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:))
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i023(:)
+
+!---     i024: H+ + O -> O+ + H
+
+!        UMIST
+
+         i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:))
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i024(:)
+
+!---     i025: CO+ + H2 -> HCO2+ + H
+
+!        UMIST
+
+         i025(:) = 9.5e-10
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i025(:)
+
+!---     i026: HCO2+ + e -> H + CO2
+
+!        UMIST
+
+         i026(:) = 1.75e-8*((300./t_elect(:))**0.5)
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i026(:)
+
+!---     i027+i028: HCO2+ + e -> H + O + CO
+
+!        UMIST
+         !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
+         !i028: 0.5 (HCO2+ + e-) -> O + CO
+
+         i027(:) = 2.38e-7*((300./t_elect(:))**0.5)
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i027(:)
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i027(:)
+
+!---     i029: HCO2+ + e -> H + CO2
+
+!        UMIST
+
+         i029(:) = 9.45e-8*((300./t_elect(:))**0.5)
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i029(:)
+
+      end if   !ionchem
 
 !----------------------------------------------------------------------
@@ -1413,5 +1438,7 @@
 !======================================================================
 
- subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
+ subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer,            &
+                        nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
+                        v_phot, v_3, v_4)
 
 !======================================================================
@@ -1420,7 +1447,4 @@
 
 use types_asis
-use photolysis_mod, only : nb_phot_max,       &
-                           nb_reaction_3_max, &
-                           nb_reaction_4_max
 
 implicit none
@@ -1431,4 +1455,10 @@
 integer             :: nesp    ! number of species in the chemistry
 integer, intent(in) :: nlayer  ! number of atmospheric layers
+integer, intent(in) :: nb_reaction_3_max 
+                               ! number of quadratic reactions
+integer, intent(in) :: nb_reaction_4_max
+                               ! number of bimolecular reactions
+integer, intent(in) :: nb_phot_max
+                               ! number of processes treated numerically as photodissociations
 
 real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
@@ -1531,8 +1561,9 @@
 !================================================================
 
- subroutine indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-                   i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-                   i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
-                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
+ subroutine indice(nb_reaction_3_max, nb_reaction_4_max,                &
+                   nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
+                   i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
+                   i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,            &
+                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
                    i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
 
@@ -1550,7 +1581,4 @@
 
 use types_asis
-use photolysis_mod, only : nb_phot_max,       &
-                           nb_reaction_3_max, &
-                           nb_reaction_4_max
 
 implicit none
@@ -1563,4 +1591,11 @@
            i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
            i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
+integer, intent(in) :: nb_reaction_3_max
+                       ! number of quadratic reactions
+integer, intent(in) :: nb_reaction_4_max
+                       ! number of bimolecular reactions
+integer, intent(in) :: nb_phot_max
+                       ! number of processes treated numerically as photodissociations
+logical, intent(in) :: ionchem
 
 ! local
@@ -1683,698 +1718,703 @@
 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n) 
 
+!Only if ion chemistry included
+if (ionchem) then
+
 !===========================================================
 !      CO2 + hv -> CO2+ + e-
 !===========================================================
 
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
+
+!===========================================================
+!      CO2 + hv -> O+ + CO + e-
+!===========================================================
+!We divide this reaction in two
+
+!0.5 CO2 + hv -> CO
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
+
+!0.5 CO2 + hv -> O+ + e-
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
+
+!===========================================================
+!      CO2 + hv -> CO+ + O + e-
+!===========================================================
+!We divide this reaction in two
+
+!0.5 CO2 + hv -> O
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
+
+!0.5 CO2 + hv -> CO+ + e-
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
+
+!===========================================================
+!      CO2 + hv -> C+ + O2 + e-
+!===========================================================
+!We divide this reaction in two
+
+!0.5 CO2 + hv -> O2
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
+
+!0.5 CO2 + hv -> C+ + e-
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
+
+!===========================================================
+!      O2 + hv -> O2+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
+
+!===========================================================
+!      O + hv -> O+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
+
+!===========================================================
+!      NO + hv -> NO+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
+
+!===========================================================
+!      CO + hv -> CO+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
+
+!===========================================================
+!      CO + hv -> C+ + O + e-
+!===========================================================
+!We divide this reaction in two
+
+!0.5 CO + hv -> O
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
+
+!0.5 CO + hv -> C+ + e-
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
+
+!===========================================================
+!      N2 + hv -> N2+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
+
+!===========================================================
+!      N2 + hv -> N+ + N + e-
+!===========================================================
+!We divide this reaction in two
+
+!0.5 N2 + hv -> N
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
+
+!0.5 N2 + hv -> N+ + e-
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
+
+!===========================================================
+!      N + hv -> N+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
+
+!===========================================================
+!      H + hv -> H+ + e-
+!===========================================================
+
+   nb_phot = nb_phot + 1
+
+   indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
+
+end if   !ionchem
+
+!===========================================================
+!      a001 : O + O2 + CO2 -> O3 + CO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy) 
+
+!===========================================================
+!      a002 : O + O + CO2 -> O2 + CO2
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      a003 : O + O3 -> O2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      b001 : O(1D) + CO2 -> O + CO2
+!===========================================================
+
 nb_phot = nb_phot + 1
 
-indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
-
-!===========================================================
-!      CO2 + hv -> O+ + CO + e-
+indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
+
+!===========================================================
+!      b002 : O(1D) + H2O -> OH + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy) 
+
+!===========================================================
+!      b003 : O(1D) + H2 -> OH + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h) 
+
+!===========================================================
+!      b004 : O(1D) + O2 -> O + O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
+
+!===========================================================
+!      b005 : O(1D) + O3 -> O2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      b006 : O(1D) + O3 -> O2 + O + O
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 
+
+!===========================================================
+!      c001 : O + HO2 -> OH + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2) 
+
+!===========================================================
+!      c002 : O + OH -> O2 + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h) 
+
+!===========================================================
+!      c003 : H + O3 -> OH + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2) 
+
+!===========================================================
+!      c004 : H + HO2 -> OH + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy) 
+
+!===========================================================
+!      c005 : H + HO2 -> H2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2) 
+
+!===========================================================
+!      c006 : H + HO2 -> H2O + O
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o) 
+
+!===========================================================
+!      c007 : OH + HO2 -> H2O + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2) 
+
+!===========================================================
+!      c008 : HO2 + HO2 -> H2O2 + O2
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
+
+!===========================================================
+!      c009 : OH + H2O2 -> H2O + HO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2) 
+
+!===========================================================
+!      c010 : OH + H2 -> H2O + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h) 
+
+!===========================================================
+!      c011 : H + O2 + CO2 -> HO2 + CO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy) 
+
+!===========================================================
+!      c012 : O + H2O2 -> OH + HO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2) 
+
+!===========================================================
+!      c013 : OH + OH -> H2O + O
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o) 
+
+!===========================================================
+!      c014 : OH + O3 -> HO2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2) 
+
+!===========================================================
+!      c015 : HO2 + O3 -> OH + O2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2) 
+
+!===========================================================
+!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
+
+!===========================================================
+!      c017 : OH + OH + CO2 -> H2O2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy) 
+
+!===========================================================
+!      c018 : H + H + CO2 -> H2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 
+
+!===========================================================
+!      d001 : NO2 + O -> NO + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2) 
+
+!===========================================================
+!      d002 : NO + O3 -> NO2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2) 
+
+!===========================================================
+!      d003 : NO + HO2 -> NO2 + OH 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh) 
+
+!===========================================================
+!      d004 : N + NO -> N2 + O 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o) 
+
+!===========================================================
+!      d005 : N + O2 -> NO + O 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o) 
+
+!===========================================================
+!      d006 : NO2 + H -> NO + OH 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh) 
+
+!===========================================================
+!      d007 : N + O -> NO
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy) 
+
+!===========================================================
+!      d008 : N + HO2 -> NO + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh) 
+
+!===========================================================
+!      d009 : N + OH -> NO + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h) 
+
+!===========================================================
+!      d010 : N(2D) + O -> N + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy) 
+
+!===========================================================
+!      d011 : N(2D) + N2 -> N + N2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy) 
+
+!===========================================================
+!      d012 : N(2D) + CO2 -> NO + CO
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co) 
+
+!===========================================================
+!      e001 : CO + OH -> CO2 + H 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h) 
+
+!===========================================================
+!      e002 : CO + O + M -> CO2 + M 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 
+
+!Only if ion chemistry
+if (ionchem) then
+
+!===========================================================
+!      i001 : CO2+ + O2 -> O2+ + CO2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
+
+!===========================================================
+!      i002 : CO2+ + O -> O+ + CO2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
+
+!===========================================================
+!      i003 : CO2+ + O -> O2+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
+
+!===========================================================
+!      i004 : O2+ + e- -> O + O 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
+
+!===========================================================
+!      i005 : O+ + CO2 -> O2+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
+
+!===========================================================
+!      i006 : CO2+ + e -> CO + O 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
+
+!===========================================================
+!      i007 : CO2+ + NO -> NO+ + CO2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
+
+!===========================================================
+!      i008 : O2+ + NO -> NO+ + O2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
+
+!===========================================================
+!      i009 : O2+ + N2 -> NO+ + NO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
+
+!===========================================================
+!      i010 : O2+ + N -> NO+ + O 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
+
+!===========================================================
+!      i011 : O+ + N2 -> NO+ + N 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
+
+!===========================================================
+!      i012 : NO+ + e -> N + O 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
+
+!===========================================================
+!      i013 : CO+ + CO2 -> CO2+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
+
+!===========================================================
+!      i014 : CO+ + O -> O+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
+
+!===========================================================
+!      i015 : C+ + CO2 -> CO+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
+
+!===========================================================
+!      i016 : N2+ + CO2 -> CO2+ + N2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
+
+!===========================================================
+!      i017 : N2+ + O -> NO+ + N 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
+
+!===========================================================
+!      i018 : N2+ + CO -> CO+ + N2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
+
+!===========================================================
+!      i019 : N2+ + e -> N + N 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
+
+!===========================================================
+!      i020 : N2+ + O -> O+ + N2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
+
+!===========================================================
+!      i021 : N+ + CO2 -> CO2+ + N 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
+
+!===========================================================
+!      i022 : CO+ + H -> H+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
+
+!===========================================================
+!      i023 : O+ + H -> H+ + O 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
+
+!===========================================================
+!      i024 : H+ + O -> O+ + H 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
+
+!===========================================================
+!      i025 : CO2+ + H2 -> HCO2+ + H 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
+
+!===========================================================
+!      i026 : HCO2+ + e -> H + CO2 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
+
+!===========================================================
+!      i027 : HCO2+ + e -> H + O + CO 
 !===========================================================
 !We divide this reaction in two
 
-!0.5 CO2 + hv -> CO
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
-
-!0.5 CO2 + hv -> O+ + e-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
-
-
-!===========================================================
-!      CO2 + hv -> CO+ + O + e-
-!===========================================================
-!We divide this reaction in two
-
-!0.5 CO2 + hv -> O
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
-
-!0.5 CO2 + hv -> CO+ + e-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
-
-
-!===========================================================
-!      CO2 + hv -> C+ + O2 + e-
-!===========================================================
-!We divide this reaction in two
-
-!0.5 CO2 + hv -> O2
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
-
-!0.5 CO2 + hv -> C+ + e-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
-
-!===========================================================
-!      O2 + hv -> O2+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
-
-!===========================================================
-!      O + hv -> O+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
-
-!===========================================================
-!      NO + hv -> NO+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
-
-!===========================================================
-!      CO + hv -> CO+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
-
-!===========================================================
-!      CO + hv -> C+ + O + e-
-!===========================================================
-!We divide this reaction in two
-
-!0.5 CO + hv -> O
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
-
-!0.5 CO + hv -> C+ + e-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
-
-
-!===========================================================
-!      N2 + hv -> N2+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
-
-!===========================================================
-!      N2 + hv -> N+ + N + e-
-!===========================================================
-!We divide this reaction in two
-
-!0.5 N2 + hv -> N
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
-
-!0.5 N2 + hv -> N+ + e-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
-
-!===========================================================
-!      N + hv -> N+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
-
-!===========================================================
-!      H + hv -> H+ + e-
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
-
-!===========================================================
-!      a001 : O + O2 + CO2 -> O3 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy) 
-
-!===========================================================
-!      a002 : O + O + CO2 -> O2 + CO2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      a003 : O + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b001 : O(1D) + CO2 -> O + CO2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b002 : O(1D) + H2O -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      b003 : O(1D) + H2 -> OH + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h) 
-
-!===========================================================
-!      b004 : O(1D) + O2 -> O + O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b005 : O(1D) + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b006 : O(1D) + O3 -> O2 + O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 
-
-!===========================================================
-!      c001 : O + HO2 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c002 : O + OH -> O2 + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h) 
-
-!===========================================================
-!      c003 : H + O3 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c004 : H + HO2 -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      c005 : H + HO2 -> H2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2) 
-
-!===========================================================
-!      c006 : H + HO2 -> H2O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c007 : OH + HO2 -> H2O + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2) 
-
-!===========================================================
-!      c008 : HO2 + HO2 -> H2O2 + O2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c009 : OH + H2O2 -> H2O + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2) 
-
-!===========================================================
-!      c010 : OH + H2 -> H2O + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h) 
-
-!===========================================================
-!      c011 : H + O2 + CO2 -> HO2 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy) 
-
-!===========================================================
-!      c012 : O + H2O2 -> OH + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2) 
-
-!===========================================================
-!      c013 : OH + OH -> H2O + O
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c014 : OH + O3 -> HO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2) 
-
-!===========================================================
-!      c015 : HO2 + O3 -> OH + O2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2) 
-
-!===========================================================
-!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c017 : OH + OH + CO2 -> H2O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy) 
-
-!===========================================================
-!      c018 : H + H + CO2 -> H2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 
-
-!===========================================================
-!      d001 : NO2 + O -> NO + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2) 
-
-!===========================================================
-!      d002 : NO + O3 -> NO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2) 
-
-!===========================================================
-!      d003 : NO + HO2 -> NO2 + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh) 
-
-!===========================================================
-!      d004 : N + NO -> N2 + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o) 
-
-!===========================================================
-!      d005 : N + O2 -> NO + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o) 
-
-!===========================================================
-!      d006 : NO2 + H -> NO + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d007 : N + O -> NO
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy) 
-
-!===========================================================
-!      d008 : N + HO2 -> NO + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d009 : N + OH -> NO + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h) 
-
-!===========================================================
-!      d010 : N(2D) + O -> N + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy) 
-
-!===========================================================
-!      d011 : N(2D) + N2 -> N + N2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy) 
-
-!===========================================================
-!      d012 : N(2D) + CO2 -> NO + CO
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co) 
-
-!===========================================================
-!      e001 : CO + OH -> CO2 + H 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h) 
-
-!===========================================================
-!      e002 : CO + O + M -> CO2 + M 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 
-
-!===========================================================
-!      i001 : CO2+ + O2 -> O2+ + CO2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
-
-!===========================================================
-!      i002 : CO2+ + O -> O+ + CO2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
-
-!===========================================================
-!      i003 : CO2+ + O -> O2+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
-
-!===========================================================
-!      i004 : O2+ + e- -> O + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
-
-!===========================================================
-!      i005 : O+ + CO2 -> O2+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
-
-!===========================================================
-!      i006 : CO2+ + e -> CO + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
-
-!===========================================================
-!      i007 : CO2+ + NO -> NO+ + CO2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
-
-!===========================================================
-!      i008 : O2+ + NO -> NO+ + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
-
-!===========================================================
-!      i009 : O2+ + N2 -> NO+ + NO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
-
-!===========================================================
-!      i010 : O2+ + N -> NO+ + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
-
-!===========================================================
-!      i011 : O+ + N2 -> NO+ + N 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
-
-!===========================================================
-!      i012 : NO+ + e -> N + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
-
-!===========================================================
-!      i013 : CO+ + CO2 -> CO2+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
-
-!===========================================================
-!      i014 : CO+ + O -> O+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
-
-!===========================================================
-!      i015 : C+ + CO2 -> CO+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
-
-!===========================================================
-!      i016 : N2+ + CO2 -> CO2+ + N2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
-
-!===========================================================
-!      i017 : N2+ + O -> NO+ + N 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
-
-!===========================================================
-!      i018 : N2+ + CO -> CO+ + N2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
-
-!===========================================================
-!      i019 : N2+ + e -> N + N 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
-
-!===========================================================
-!      i020 : N2+ + O -> O+ + N2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
-
-!===========================================================
-!      i021 : N+ + CO2 -> CO2+ + N 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
-
-!===========================================================
-!      i022 : CO+ + H -> H+ + CO 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
-
-!===========================================================
-!      i023 : O+ + H -> H+ + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
-
-!===========================================================
-!      i024 : H+ + O -> O+ + H 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
-
-!===========================================================
-!      i025 : CO2+ + H2 -> HCO2+ + H 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
-
-!===========================================================
-!      i026 : HCO2+ + e -> H + CO2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
-
-!===========================================================
-!      i027 : HCO2+ + e -> H + O + CO 
-!===========================================================
-!We divide this reaction in two
-
 !0.5HCO2+ + 0.5e -> H
 
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
 
 !0.5 HCO2+ + 0.5 e -> O + CO
 
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
 
 !===========================================================
@@ -2382,7 +2422,9 @@
 !===========================================================
 
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
+   nb_reaction_4 = nb_reaction_4 + 1
+
+   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
+
+end if    !ionchem
 
 !===========================================================
@@ -2455,5 +2497,5 @@
 !*****************************************************************
 
-      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
+      subroutine gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,&
                            i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
                            i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
@@ -2463,5 +2505,4 @@
                            i_hplus, i_hco2plus, i_elec, dens, rm, c) 
         
-
 !*****************************************************************
 
@@ -2485,4 +2526,5 @@
       integer, intent(in) :: nlayer ! number of atmospheric layers
       integer, intent(in) :: nq     ! number of tracers in the gcm
+      logical, intent(in) :: ionchem
       integer :: nesp               ! number of species in the chemistry
       integer :: lswitch            ! interface level between chemistries
@@ -2510,5 +2552,4 @@
       integer      :: l, iesp
       logical,save :: firstcall = .true.
-      
       
 !     first call initializations
@@ -2586,48 +2627,50 @@
             stop
          end if
-         if (igcm_co2plus == 0) then
-            write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
-            stop
-         end if
-         if (igcm_oplus == 0) then
-            write(*,*) "gcmtochim: Error; no O+ tracer !!!"
-            stop
-         end if
-         if (igcm_o2plus == 0) then
-            write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
-            stop
-         end if
-         if (igcm_noplus == 0) then
-            write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
-            stop
-         endif
-         if (igcm_coplus == 0) then
-            write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
-            stop
-         endif
-         if (igcm_cplus == 0) then
-            write(*,*) "gcmtochim: Error; no C+ tracer !!!"
-            stop
-         endif
-         if (igcm_n2plus == 0) then
-            write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
-            stop
-         endif
-         if (igcm_nplus == 0) then
-            write(*,*) "gcmtochim: Error; no N+ tracer !!!"
-            stop
-         endif
-         if (igcm_hplus == 0) then
-            write(*,*) "gcmtochim: Error; no H+ tracer !!!"
-            stop
-         endif
-         if (igcm_hco2plus == 0) then
-            write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
-            stop
-         endif
-         if (igcm_elec == 0) then
-            write(*,*) "gcmtochim: Error; no e- tracer !!!"
-            stop
-         end if
+         if (ionchem) then
+            if (igcm_co2plus == 0) then
+               write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
+               stop
+            end if
+            if (igcm_oplus == 0) then
+               write(*,*) "gcmtochim: Error; no O+ tracer !!!"
+               stop
+            end if
+            if (igcm_o2plus == 0) then
+               write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
+               stop
+            end if
+            if (igcm_noplus == 0) then
+               write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
+               stop
+            endif
+            if (igcm_coplus == 0) then
+               write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
+               stop
+            endif
+            if (igcm_cplus == 0) then
+               write(*,*) "gcmtochim: Error; no C+ tracer !!!"
+               stop
+            endif
+            if (igcm_n2plus == 0) then
+               write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
+               stop
+            endif
+            if (igcm_nplus == 0) then
+               write(*,*) "gcmtochim: Error; no N+ tracer !!!"
+               stop
+            endif
+            if (igcm_hplus == 0) then
+               write(*,*) "gcmtochim: Error; no H+ tracer !!!"
+               stop
+            endif
+            if (igcm_hco2plus == 0) then
+               write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
+               stop
+            endif
+            if (igcm_elec == 0) then
+               write(*,*) "gcmtochim: Error; no e- tracer !!!"
+               stop
+            end if
+         end if  ! ionchem
          firstcall = .false.
       end if ! of if (firstcall)
@@ -2655,16 +2698,21 @@
          rm(l,i_no2)  = zycol(l, igcm_no2)
          rm(l,i_n2)   = zycol(l, igcm_n2)
-         rm(l,i_co2plus) = zycol(l, igcm_co2plus)
-         rm(l,i_oplus) = zycol(l, igcm_oplus)
-         rm(l,i_o2plus) = zycol(l, igcm_o2plus)
-         rm(l,i_noplus) = zycol(l,igcm_noplus)
-         rm(l,i_coplus) = zycol(l,igcm_coplus)
-         rm(l,i_cplus)  = zycol(l,igcm_cplus)
-         rm(l,i_n2plus) = zycol(l,igcm_n2plus)
-         rm(l,i_nplus)  = zycol(l,igcm_nplus)
-         rm(l,i_hplus)  = zycol(l,igcm_hplus)
-         rm(l,i_hco2plus)=zycol(l,igcm_hco2plus)
-         rm(l,i_elec) = zycol(l, igcm_elec)
       end do 
+
+      if (ionchem) then
+         do l = 1,nlayer
+            rm(l,i_co2plus)  = zycol(l, igcm_co2plus)
+            rm(l,i_oplus)    = zycol(l, igcm_oplus)
+            rm(l,i_o2plus)   = zycol(l, igcm_o2plus)
+            rm(l,i_noplus)   = zycol(l, igcm_noplus)
+            rm(l,i_coplus)   = zycol(l, igcm_coplus)
+            rm(l,i_cplus)    = zycol(l, igcm_cplus)
+            rm(l,i_n2plus)   = zycol(l, igcm_n2plus)
+            rm(l,i_nplus)    = zycol(l, igcm_nplus)
+            rm(l,i_hplus)    = zycol(l, igcm_hplus)
+            rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
+            rm(l,i_elec)     = zycol(l, igcm_elec)
+         end do 
+      end if
 
       where (rm(:,:) < 1.e-30)
@@ -2686,12 +2734,11 @@
 !*****************************************************************
  
-      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,         &
-                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
-                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
-                           i_n, i_n2d, i_no, i_no2, i_n2,            &
-                           i_co2plus, i_oplus, i_o2plus, i_noplus,   &
-                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
+      subroutine chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, &
+                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
+                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
+                           i_n, i_n2d, i_no, i_no2, i_n2,             &
+                           i_co2plus, i_oplus, i_o2plus, i_noplus,    &
+                           i_coplus, i_cplus, i_n2plus, i_nplus,      &
                            i_hplus, i_hco2plus, i_elec, dens, c) 
-        
  
 !*****************************************************************
@@ -2716,4 +2763,5 @@
       integer, intent(in) :: nlayer  ! number of atmospheric layers
       integer, intent(in) :: nq      ! number of tracers in the gcm
+      logical, intent(in) :: ionchem
       integer :: nesp                ! number of species in the chemistry
       integer :: lswitch             ! interface level between chemistries
@@ -2762,19 +2810,23 @@
          zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
          zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
-         zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
-         zycol(l, igcm_oplus)   = c(l,i_oplus)/dens(l)
-         zycol(l, igcm_o2plus)  = c(l,i_o2plus)/dens(l)
-         zycol(l, igcm_noplus)  = c(l,i_noplus)/dens(l)
-         zycol(l, igcm_coplus)  = c(l,i_coplus)/dens(l)
-         zycol(l, igcm_cplus)   = c(l,i_cplus)/dens(l)
-         zycol(l, igcm_n2plus)  = c(l,i_n2plus)/dens(l)
-         zycol(l, igcm_nplus)   = c(l,i_nplus)/dens(l)
-         zycol(l, igcm_hplus)   = c(l,i_hplus)/dens(l)
-         zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
-         zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
       end do 
 
+      if (ionchem) then
+         do l = 1,lswitch-1
+            zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
+            zycol(l, igcm_oplus)   = c(l,i_oplus)/dens(l)
+            zycol(l, igcm_o2plus)  = c(l,i_o2plus)/dens(l)
+            zycol(l, igcm_noplus)  = c(l,i_noplus)/dens(l)
+            zycol(l, igcm_coplus)  = c(l,i_coplus)/dens(l)
+            zycol(l, igcm_cplus)   = c(l,i_cplus)/dens(l)
+            zycol(l, igcm_n2plus)  = c(l,i_n2plus)/dens(l)
+            zycol(l, igcm_nplus)   = c(l,i_nplus)/dens(l)
+            zycol(l, igcm_hplus)   = c(l,i_hplus)/dens(l)
+            zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
+            zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
+         end do
+      end if 
+
       end subroutine chimtogcm
 
 end subroutine photochemistry
-
Index: trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 2164)
+++ trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 2170)
@@ -1,5 +1,5 @@
 !==========================================================================
 
-      subroutine photolysis(nlayer,                              &
+      subroutine photolysis(nlayer, nb_phot_max,                  &
                             lswitch, press, temp, sza, tauref,   &
                             zmmean, dist_sol, rmco2, rmo3, v_phot)
@@ -8,5 +8,4 @@
 
       use comcstfi_h
-      use photolysis_mod, only : nb_phot_max
 
       implicit none
@@ -18,14 +17,15 @@
 !==========================================================================
        
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer :: lswitch            ! interface level between chemistries
-      real :: press(nlayer)         ! pressure (hPa)
-      real :: temp(nlayer)          ! temperature (K)
-      real :: sza                   ! solar zenith angle (deg)
-      real :: tauref                ! optical depth at 7 hpa
-      real :: zmmean(nlayer)        ! mean molecular mass (g)
-      real :: dist_sol              ! sun distance (AU) 
-      real :: rmco2(nlayer)         ! co2 mixing ratio
-      real :: rmo3(nlayer)          ! ozone mixing ratio
+      integer, intent(in) :: nlayer      ! number of atmospheric layers
+      integer, intent(in) :: nb_phot_max ! number of processes treated numerically as photodissociations
+      integer :: lswitch                 ! interface level between chemistries
+      real :: press(nlayer)              ! pressure (hPa)
+      real :: temp(nlayer)               ! temperature (K)
+      real :: sza                        ! solar zenith angle (deg)
+      real :: tauref                     ! optical depth at 7 hpa
+      real :: zmmean(nlayer)             ! mean molecular mass (g)
+      real :: dist_sol                   ! sun distance (AU) 
+      real :: rmco2(nlayer)              ! co2 mixing ratio
+      real :: rmo3(nlayer)               ! ozone mixing ratio
 
 !==========================================================================
Index: trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90	(revision 2164)
+++ trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90	(revision 2170)
@@ -4,14 +4,7 @@
 
 ! photolysis
-  logical, parameter :: jonline = .true.  ! true: on-line ! false: lookup table
+
   integer, parameter :: nphot = 13        ! number of photolysis
-  integer, parameter :: nphotion = 18     ! number of photoionizations
   integer, parameter :: nabs  = 10        ! number of absorbing gases
-
-! number of reactions in chemical solver
-
-  integer, parameter :: nb_phot_max       = nphot + nphotion + 9 ! photolysis + photoionization + quenching/heterogeneous
-  integer, parameter :: nb_reaction_3_max = 6         ! quadratic
-  integer, parameter :: nb_reaction_4_max = 60        ! bimolecular
 
 ! spectral grid
Index: trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F	(revision 2164)
+++ trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F	(revision 2170)
@@ -1,6 +1,6 @@
 !==============================================================================
 
-      subroutine photolysis_online(nlayer, alt, press, temp,
-     $           zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
+      subroutine photolysis_online(nlayer, nb_phot_max, alt, press, 
+     $           temp, zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
      $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
      $           i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
@@ -15,4 +15,5 @@
       integer, intent(in) :: nesp                                    ! total number of chemical species
       integer, intent(in) :: nlayer
+      integer, intent(in) :: nb_phot_max
       integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
      $                       i_h2, i_oh, i_ho2, i_h2o2, i_h2o, 
