Index: trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2316)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2321)
@@ -23,5 +23,5 @@
                             igcm_noplus, igcm_n2plus, igcm_hplus,         &
                             igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
-                            igcm_elec, mmol
+                            igcm_h3oplus, igcm_ohplus, igcm_elec, mmol
 
       use conc_mod, only: mmean ! mean molecular mass of the atmosphere
@@ -144,4 +144,6 @@
       integer,save :: i_hcoplus=0
       integer,save :: i_h2oplus=0
+      integer,save :: i_h3oplus=0
+      integer,save :: i_ohplus=0
       integer,save :: i_elec=0
 
@@ -577,4 +579,36 @@
             if (i_h2oplus /= 0) then
                write(*,*) "calchim: Error: H2O+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
+         i_h3oplus=igcm_h3oplus
+         if(ionchem) then
+            if (i_h3oplus == 0) then
+               write(*,*) "calchim: Error, no H3O+ tracer !!!"
+               write(*,*) "H3O+ is needed if O2+ is in traceur.def"
+               call abort_physic("calchim","missing h3oplus tracer",1)
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_h3oplus
+            end if
+         else
+            if (i_h3oplus /= 0) then
+               write(*,*) "calchim: Error: H3O+ is present, but O2+ is not!!!"
+               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
+            endif
+         endif
+         i_ohplus=igcm_ohplus
+         if(ionchem) then
+            if (i_ohplus == 0) then
+               write(*,*) "calchim: Error, no OH+ tracer !!!"
+               write(*,*) "OH+ is needed if O2+ is in traceur.def"
+               call abort_physic("calchim","missing ohplus tracer",1)
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_ohplus
+            end if
+         else
+            if (i_ohplus /= 0) then
+               write(*,*) "calchim: Error: OH+ is present, but O2+ is not!!!"
                write(*,*) "Both must be in traceur.def if ion chemistry wanted"
             endif
@@ -665,5 +699,5 @@
             ! set number of reactions, depending on ion chemistry or not
             if (ionchem) then
-               nb_reaction_4_max = 80   ! set number of bimolecular reactions
+               nb_reaction_4_max = 95   ! set number of bimolecular reactions
                nb_reaction_3_max = 6    ! set number of quadratic reactions
                nquench           = 9    ! set number of quenching + heterogeneous reactions
@@ -688,5 +722,5 @@
                            zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
                            dist_sol,zday,surfdust1d,surfice1d,            &
-                           jo3,jh2o,taucol,iter)
+                           jo3,jh2o,em_no,em_o2,taucol,iter)
 
 !        ozone photolysis, for output
@@ -711,10 +745,10 @@
                     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)
-               end do
-            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  ! photochem
 
@@ -773,10 +807,10 @@
                              ' ',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)
                call writediagfi(ngrid,'emission_o2',        &
                     'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
-            endif
+!            endif
             
             if (callstats) then
Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 2316)
+++ trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 2321)
@@ -87,6 +87,4 @@
       igcm_h2o_vap=0
       igcm_h2o_ice=0
-      igcm_hdo_vap=0
-      igcm_hdo_ice=0
       igcm_co2=0
       igcm_co=0
@@ -121,4 +119,6 @@
       igcm_hcoplus=0
       igcm_h2oplus=0
+      igcm_h3oplus=0
+      igcm_ohplus=0
       igcm_elec=0
 
@@ -278,14 +278,4 @@
            count = count + 1
         end if
-        if (noms(iq) == "hdo_vap") then
-           igcm_hdo_vap = iq
-           mmol(igcm_hdo_vap) = 18.
-           count = count + 1
-        end if
-        if (noms(iq) == "hdo_ice") then
-           igcm_hdo_ice = iq
-           mmol(igcm_hdo_ice) = 18.
-           count = count + 1
-        end if
         if (noms(iq).eq."he") then
           igcm_he=iq
@@ -354,4 +344,14 @@
            igcm_h2oplus = iq
            mmol(igcm_h2oplus) = 18.
+           count = count + 1
+        end if
+        if (noms(iq) == "h3oplus") then
+           igcm_h3oplus = iq
+           mmol(igcm_h3oplus) = 19.
+           count = count + 1
+        end if
+        if (noms(iq) == "ohplus") then
+           igcm_ohplus = iq
+           mmol(igcm_ohplus) = 17.
            count = count + 1
         end if
@@ -573,5 +573,5 @@
                   do n = 1,nspe
                      iq = niq(n)
-                     if (iq /= igcm_h2o_vap .or. iq/= igcm_hdo_vap .or. flagh2o == 1) then
+                     if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
                         pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
                      end if
@@ -617,9 +617,10 @@
               .or. igcm_cplus == 0 .or. igcm_nplus == 0 .or. igcm_noplus == 0    & 
               .or. igcm_n2plus == 0 .or. igcm_hplus == 0 .or. igcm_hco2plus == 0 &
-              .or. igcm_hcoplus == 0 .or. igcm_h2oplus == 0 .or. igcm_elec == 0) then
+              .or. igcm_hcoplus == 0 .or. igcm_h2oplus == 0                      &
+              .or. igcm_h3oplus == 0.or. igcm_ohplus == 0 .or.igcm_elec == 0) then
             write(*,*)'inichim_newstart error:'
             write(*,*)'if co2plus is in traceur.def, all other ions must also be'
             write(*,*)'o2plus, oplus, coplus, cplus, nplus, noplus, n2plus'
-            write(*,*)'hplus, hco2plus, hcoplus, h2oplus, and elec'
+            write(*,*)'hplus, hco2plus, hcoplus, h2oplus, h3oplus, ohplus, and elec'
             write(*,*)'stop'
             call abort_physic("inichim_newstart","missing ions in tracers",1)
@@ -642,4 +643,6 @@
                   pq(i,j,l,igcm_hcoplus)  = 0.
                   pq(i,j,l,igcm_h2oplus)  = 0.
+                  pq(i,j,l,igcm_h3oplus)  = 0.
+                  pq(i,j,l,igcm_ohplus)   = 0.
                   pq(i,j,l,igcm_elec)     = 0.
                end do
@@ -661,4 +664,6 @@
          qsurf(1:ngrid,igcm_hcoplus)  = 0.
          qsurf(1:ngrid,igcm_h2oplus)  = 0.
+         qsurf(1:ngrid,igcm_h3oplus)  = 0.
+         qsurf(1:ngrid,igcm_ohplus)   = 0.
          qsurf(1:ngrid,igcm_elec)     = 0.
 
@@ -668,5 +673,6 @@
               .or. igcm_cplus /= 0 .or. igcm_nplus /= 0 .or. igcm_noplus /= 0    & 
               .or. igcm_n2plus /= 0 .or. igcm_hplus /= 0 .or. igcm_hco2plus /= 0 &
-              .or. igcm_hcoplus /= 0 .or. igcm_h2oplus /= 0 .or. igcm_elec /= 0) then
+              .or. igcm_hcoplus /= 0 .or. igcm_h2oplus /= 0                      &
+              .or. igcm_h3oplus /= 0 .or. igcm_ohplus /= 0 .or. igcm_elec /= 0) then
             write(*,*)'inichim_newstart error:'
             write(*,*)'some ions are in traceur.def, but not co2plus'
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2316)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2321)
@@ -18,5 +18,6 @@
                           alt, temp, temp_elect, dens, zmmean,                 &
                           dist_sol, zday,                                      &
-                          surfdust1d, surfice1d, jo3, jh2o,tau, iter)
+                          surfdust1d, surfice1d, jo3, jh2o,em_no,em_o2,        &
+                          tau, iter)
 
 use param_v4_h, only: jion
@@ -72,4 +73,6 @@
 real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
 real    :: jh2o(nlayer)        ! photodissociation rate h2o -> h + oh
+real    :: em_no(nlayer)       ! NO nightglow volume emission rate
+real    :: em_o2(nlayer)       ! O2 nightglow volume emission rate
 
 !===================================================================
@@ -114,5 +117,7 @@
 integer,parameter :: i_hcoplus = 28
 integer,parameter :: i_h2oplus = 29
-integer,parameter :: i_elec    = 30
+integer,parameter :: i_h3oplus = 30
+integer,parameter :: i_ohplus  = 31
+integer,parameter :: i_elec    = 32
 
 integer :: ilay
@@ -157,5 +162,5 @@
                i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
                i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
-               i_h2oplus, i_elec)
+               i_h2oplus, i_h3oplus, i_ohplus, i_elec)
    firstcall = .false.
 end if
@@ -171,5 +176,6 @@
                i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
                i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
-               i_hcoplus, i_h2oplus, i_elec, dens, rm, c) 
+               i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus,      &
+               i_elec, dens, rm, c) 
 
 !===================================================================
@@ -211,6 +217,4 @@
          v_phot(:,31)=jion(12,:,1)
       endif
-!      write(*,*)'photochemistry/205',c(:,i_co2),ig
-!      write(*,*)'photochemistry/206',v_phot(:,3),ig
       
    else ! night
@@ -367,13 +371,16 @@
            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)+                &
-           c(ilev,i_hcoplus)+c(ilev,i_h2oplus)) then
+           c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+c(ilev,i_h3oplus)+             &
+           c(ilev,i_ohplus)) 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)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)
+              c(ilev,i_hco2plus)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+         &
+              c(ilev,i_h3oplus)+c(ilev,i_ohplus)
          !      write(*,*)'photochemistry/359'
          !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
       end if
    end if
+
    cnew(:)   = 0.
 
@@ -386,4 +393,8 @@
 
 end do ! ilev
+
+!add calculation of NO and O2 nightglows
+em_no(:)=c(:,i_o)*c(:,i_n)*v_4(:,26)   !2.8e-17*(300./temp(:)))**0.5
+em_o2(:)=0.75*c(:,i_o)*c(:,i_o)*c(:,i_co2)*v_3(:,1)   !2.5*9.46e-34*exp(485./temp(:))*dens(:)
 
 !===================================================================
@@ -397,5 +408,6 @@
                i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
                i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
-               i_hcoplus, i_h2oplus, i_elec, dens, c) 
+               i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus,      &
+               i_elec, dens, c) 
 contains
 
@@ -599,5 +611,7 @@
                            i034, i035, i036, i037, i038, i039, i040,   &
                            i041, i042, i043, i044, i045, i046, i047,   &
-                           i048, i049,                                 &
+                           i048, i049, i050, i051, i052, i053, i054,   &
+                           i055, i056, i057, i058, i059, i060, i061,   &
+                           i062, i063,                                 &
                            h001, h002, h003, h004, h005
 
@@ -1521,4 +1535,120 @@
          nb_reaction_4 = nb_reaction_4 + 1
          v_4(:,nb_reaction_4) = i049(:)
+
+!---     i050: H2O+ + H2O -> H3O+ + OH
+
+         !UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973)
+
+         i050(:) = 2.1e-9*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i050(:)
+
+
+!---     i051: H2O+ + H2 -> H3O+ + H
+
+         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
+
+         i051(:) = 6.4e-10
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i051(:)
+
+!---     i052: HCO+ + H2O -> H3O+ + CO
+
+         !UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978)
+
+         i052(:) = 2.5e-9*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i052(:)
+
+!---     i053: H3O+ + e -> OH + H + H
+
+         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
+
+         i053(:) = 3.05e-7*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i053(:)
+
+!---     i054: H3O + + e -> H2O + H
+
+         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
+         
+         i054(:) = 7.09e-8*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i054(:)
+
+!---     i055: H3O+ + e -> OH + H2
+
+         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
+
+         i055(:) = 5.37e-8*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i055(:)
+
+!---     i056: H3O+ + e -> O + H2 + H
+
+         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
+
+         i056(:) = 5.6e-9*((300./t_elect(:))**0.5)
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i056(:)
+
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i056(:)
+
+!---     i057: O+ + H2 -> OH+ + H 
+
+         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
+
+         i057(:) = 1.7e-9
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i057(:)
+
+!---     i058: OH+ + O -> O2+ + H
+
+         !UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1
+
+         i058(:) = 7.1e-10
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i058(:)
+
+!---     i059: OH+ + CO2 -> HCO2+ + O
+
+         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
+
+         i059(:) = 1.44e-9
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i059(:)
+
+!---     i060: OH+ + CO -> HCO+ + O
+
+         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
+
+         i060(:) = 1.05e-9
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i060(:)
+
+!---     i061: OH+ + NO -> NO+ + OH (tasa de reacciÃ³n UMIST 3.59e-10)
+
+         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
+
+         i061(:) = 3.59e-10
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i061(:)
+
+!---     i062: OH+ + H2 -> H2O+ + H (tasa de reacciÃ³n UMIST 1.01e-9,
+
+         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
+
+         i062(:) = 1.01e-9
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i062(:)
+
+!---     i063: OH+ + O2 -> O2+ + OH (tasa de reacciÃ³n UMIST 5.9e-10
+
+         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
+
+         i063(:) = 5.9e-10
+         nb_reaction_4 = nb_reaction_4 + 1
+         v_4(:,nb_reaction_4) = i063(:)
 
       end if   !ionchem
@@ -1734,5 +1864,5 @@
                    i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
                    i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
-                   i_h2oplus, i_elec)
+                   i_h2oplus, i_h3oplus, i_ohplus, i_elec)
 
 !================================================================
@@ -1759,5 +1889,5 @@
            i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
            i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
-           i_hcoplus, i_h2oplus, i_elec
+           i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec
 integer, intent(in) :: nb_reaction_3_max
                        ! number of quadratic reactions
@@ -2744,4 +2874,110 @@
    indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
 
+!===========================================================
+!      i050 : H2O+ + H2O -> H3O+ + OH 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh)
+
+!===========================================================
+!      i051 : H2O+ + H2 -> H3O+ + H 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h)
+
+!===========================================================
+!      i052 : HCO+ + H2O -> H3O+ + CO 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co)
+
+!===========================================================
+!      i053: H3O+ + e -> OH + H + H
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 2.0, i_h)
+
+!===========================================================
+!      i054: H3O+ + e -> H2O + H
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_h2o, 1.0, i_h)
+
+!===========================================================
+!      i055: H3O+ + e -> HO + H2
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2)
+
+!===========================================================
+!      i056: H3O+ + e -> O + H2 + H
+!===========================================================
+!We divide this reaction in two
+
+!0.5H3O+ + 0.5e -> O
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy)
+
+!0.5H3O+ + 0.5e -> H2 + H
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h)
+
+!===========================================================
+!      i057: O+ + H2 -> OH+ + H
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h)
+
+!===========================================================
+!      i058: OH+ + O -> O2+ + H
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h)
+
+!===========================================================
+!      i059: OH+ + CO2 -> HCO2+ + O
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co2, 1.0, i_hco2plus, 1.0, i_o)
+
+!===========================================================
+!      i060: OH+ + CO -> HCO+ + O
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o)
+
+!===========================================================
+!      i061: OH+ + NO -> NO+ + OH 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh)
+
+!===========================================================
+!      i062: OH+ + H2 -> H2O+ + H 
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h)
+
+!===========================================================
+!      i063: OH+ + O2 -> O2+ + OH
+!===========================================================
+
+   nb_reaction_4 = nb_reaction_4 + 1
+   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh)
+
 end if    !ionchem
 
@@ -2822,5 +3058,5 @@
                            i_coplus, i_cplus, i_n2plus, i_nplus,     &
                            i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,&
-                           i_elec, dens, rm, c) 
+                           i_h3oplus, i_ohplus, i_elec, dens, rm, c) 
         
 !*****************************************************************
@@ -2834,5 +3070,5 @@
      &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
      &                      igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,   &
-     &                      igcm_elec
+     &                      igcm_h3oplus, igcm_ohplus, igcm_elec
 
       implicit none
@@ -2855,5 +3091,5 @@
                  i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus,   &
                  i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus,    &
-                 i_hcoplus, i_h2oplus, i_elec
+                 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec
 
       real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
@@ -2997,4 +3233,12 @@
                call abort_physic("gcmtochim","missing h2oplus tracer",1)
             endif
+            if (igcm_h3oplus == 0) then
+               write(*,*) "gcmtochim: Error; no H3O+ tracer !!!"
+               call abort_physic("gcmtochim","missing h3oplus tracer",1)
+            endif
+            if (igcm_ohplus == 0) then
+               write(*,*) "gcmtochim: Error; no OH+ tracer !!!"
+               call abort_physic("gcmtochim","missing ohplus tracer",1)
+            endif
             if (igcm_elec == 0) then
                write(*,*) "gcmtochim: Error; no e- tracer !!!"
@@ -3043,4 +3287,6 @@
             rm(l,i_hcoplus)  = zycol(l, igcm_hcoplus)
             rm(l,i_h2oplus)  = zycol(l, igcm_h2oplus)
+            rm(l,i_h3oplus)  = zycol(l, igcm_h3oplus)
+            rm(l,i_ohplus)   = zycol(l, igcm_ohplus)
             rm(l,i_elec)     = zycol(l, igcm_elec)
          end do 
@@ -3072,5 +3318,5 @@
                            i_coplus, i_cplus, i_n2plus, i_nplus,      &
                            i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
-                           i_elec, dens, c) 
+                           i_h3oplus, i_ohplus, i_elec, dens, c) 
  
 !*****************************************************************
@@ -3084,5 +3330,5 @@
                             igcm_n2plus, igcm_nplus, igcm_hplus,          &
                             igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
-                            igcm_elec
+                            igcm_h3oplus, igcm_ohplus, igcm_elec
 
       implicit none
@@ -3104,5 +3350,6 @@
                  i_co2plus, i_oplus, i_o2plus, i_noplus,         &
                  i_coplus, i_cplus, i_n2plus, i_nplus,           &
-                 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, i_elec
+                 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,      &
+                 i_h3oplus, i_ohplus, i_elec
 
       real :: dens(nlayer)     ! total number density (molecule.cm-3) 
@@ -3159,4 +3406,6 @@
             zycol(l, igcm_hcoplus) = c(l,i_hcoplus)/dens(l)
             zycol(l, igcm_h2oplus) = c(l,i_h2oplus)/dens(l)
+            zycol(l, igcm_h3oplus) = c(l,i_h3oplus)/dens(l)
+            zycol(l, igcm_ohplus)  = c(l,i_ohplus)/dens(l)
             zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
          end do
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2316)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2321)
@@ -2,5 +2,5 @@
 
        use tracer_mod
-       USE comcstfi_h
+       use comcstfi_h, only: pi
        IMPLICIT NONE
 c=======================================================================
@@ -108,4 +108,6 @@
       igcm_hcoplus=0
       igcm_h2oplus=0
+      igcm_h3oplus=0
+      igcm_ohplus=0
       igcm_elec=0
 
@@ -335,4 +337,14 @@
           igcm_h2oplus=iq
           mmol(igcm_h2oplus)=18.
+          count=count+1
+        endif
+        if (noms(iq).eq."h3oplus") then
+          igcm_h3oplus=iq
+          mmol(igcm_h3oplus)=19.
+          count=count+1
+        endif
+        if (noms(iq).eq."ohplus") then
+          igcm_ohplus=iq
+          mmol(igcm_ohplus)=17.
           count=count+1
         endif
Index: trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2316)
+++ trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2321)
@@ -90,4 +90,6 @@
       integer,save :: igcm_hcoplus
       integer,save :: igcm_h2oplus
+      integer,save :: igcm_h3oplus
+      integer,save :: igcm_ohplus
       integer,save :: igcm_elec
       ! other tracers
