Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2157)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2158)
@@ -1,3 +1,3 @@
-      subroutine calchim(ngrid,nlayer,nq,                                   &
+      subroutine calchim(ngrid,nlayer,nq,                           &
                          ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
                          zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
@@ -17,4 +17,5 @@
       use comcstfi_h
       use photolysis_mod
+      use iono_h, only: temp_elect
 
       implicit none
@@ -161,4 +162,5 @@
       real :: zdens(nlayer)        ! Density  (cm-3)
       real :: ztemp(nlayer)        ! Temperature (K)
+      real :: ztelec(nlayer)       ! Electronic temperature (K)
       real :: zlocal(nlayer)       ! Altitude (km)
       real :: zycol(nlayer,nq)     ! Composition (mole fractions)
@@ -186,7 +188,16 @@
                                     !  within one physical timestep
 
+      logical :: unichim            ! only one, unified chemical scheme at all 
+                                    ! layers (default), or upper atmospheric 
+                                    ! scheme in the thermosphere?
+
 !=======================================================================
 !     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
@@ -201,9 +212,16 @@
             end if
          end if
+
+         if(.not.unichim) then
+            !Read reaction rates from external file if the upper atmospheric
+            !chemistry is called
+            call chemthermos_readini
+         endif
+
          ! find index of chemical tracers to use
          allocate(niq(nq))
          ! Listed here are all tracers that can go into photochemistry
          nbq = 0 ! to count number of tracers
-         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
+         ! Species ALWAYS present if photochem=.T.
          i_co2 = igcm_co2
          if (i_co2 == 0) then
@@ -358,241 +376,93 @@
             niq(nbq) = i_h2o
          end if
-         !Check tracers needed for thermospheric chemistry
-         if(thermochem) then
-            chemthermod=0  !Default: C/O/H chemistry
-            !Nitrogen chemistry
-            !NO is used to determine if N chemistry is wanted
-            !chemthermod=2 -> N chemistry
-            if (i_no == 0) then
-               write(*,*) "calchim: no NO tracer"
-               write(*,*) "C/O/H themosp chemistry only "
-            else
-               chemthermod=2
-               write(*,*) "calchim: NO in traceur.def"
-               write(*,*) "Nitrogen chemistry included"
-            end if
-            ! N
-            if(chemthermod == 2) then
-               if (i_n == 0) then
-                  write(*,*) "calchim: Error; no N tracer !!!"
-                  write(*,*) "N is needed if NO is in traceur.def"
-                  stop
-               end if
-            ! NO2
-               if (i_no2 == 0) then
-                  write(*,*) "calchim: Error; no NO2 tracer !!!"
-                  write(*,*) "NO2 is needed if NO is in traceur.def"
-                  stop
-               end if
-            ! N(2D)
-               if (i_n2d == 0) then
-                  write(*,*) "calchim: Error; no N2D !!!"
-                  write(*,*) "N2D is needed if NO is in traceur.def"
-                  stop
-               end if
-            endif    !Of if(chemthermod == 2)
-            ! Ions
-            ! O2+ is used to determine if ion chemistry is needed
-            ! chemthermod=3 -> ion chemistry
-            i_o2plus = igcm_o2plus
-            if(chemthermod == 2) then
-               if (i_o2plus == 0) then
-                  write(*,*) "calchim: no O2+ tracer; no ion chemistry"
-               else
-                  nbq = nbq + 1
-                  niq(nbq) = i_o2plus
-                  chemthermod = 3
-                  write(*,*) "calchim: O2+ in traceur.def"
-                  write(*,*) "Ion chemistry included"
-               end if
-            else
-               if (i_o2plus /= 0) then
-                  write(*,*) "calchim: O2+ is present, but NO is not!!!"
-                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if(chemthermod == 2)
-            ! CO2+
-            i_co2plus = igcm_co2plus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif    !Of if(chemthermod == 3)
-            ! O+
-            i_oplus = igcm_oplus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-            ! CO+
-            i_coplus = igcm_coplus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   ! Of if (chemthermod == 3)
-            ! C+
-            i_cplus = igcm_cplus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   ! Of if (chemthermod == 3)
-            ! N+
-            i_nplus = igcm_nplus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-            ! NO+
-            i_noplus = igcm_noplus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-            ! N2+
-            i_n2plus = igcm_n2plus
-            if (chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-            !H+
-            i_hplus = igcm_hplus
-            if (chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-            ! HCO2+
-            i_hco2plus = igcm_hco2plus
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif    !Of if(chemthermod == 3)
-            !e-
-            i_elec = igcm_elec
-            if(chemthermod == 3) 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 ionosphere wanted"
-                  stop
-               endif
-            endif   !Of if (chemthermod == 3)
-         endif      !Of thermochem
-
+         i_o2plus = igcm_o2plus
+         if (i_o2plus == 0) then
+            write(*,*) "calchim: Error, no O2+ tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_o2plus
+         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
+         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
+         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
+         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
+         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
+         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
+         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
+         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
+         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
+         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
+         
          write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
                
@@ -634,4 +504,6 @@
             zlocal(l) = zzlay(ig,l)/1000.
             zmmean(l) = mmean(ig,l)
+            ztelec(l) = temp_elect(zlocal(l),ztemp(l),1)
+            !Electronic temperature. Index 1 -> Viking; Index 2-> MAVEN
 
 !           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
@@ -640,19 +512,13 @@
             surfice1d(l)  = surfice(ig,l)*1.e-2
 
-!           search for switch index between regions
-
-            if (photochem .and. thermochem) then
+!           search for switch index between regions  
+            if(unichim) then
+               lswitch=nlayer+1
+            else
                if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
                   lswitch = l
                   foundswitch = 1
                end if
-            end if
-            if (.not. photochem) then
-               lswitch = 22
-            end if
-            if (.not. thermochem) then
-               lswitch = min(50,nlayer+1)
-            end if
-
+            endif
          end do ! of do l=1,nlayer
 
@@ -664,11 +530,12 @@
 !=======================================================================
 
-!        chemistry in lower atmosphere
-
+!        chemistry: only one scheme at all layers
+         
          if (photochem) then
             call photochemistry(nlayer,nq,                            &
                            ig,lswitch,zycol,szacol,ptimestep,         &
-                           zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
-                           zday,surfdust1d,surfice1d,jo3,jh2o,taucol,iter)
+                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,   &
+                           dist_sol,zday,surfdust1d,surfice1d,        &
+                           jo3,jh2o,taucol,iter)
 
 !        ozone photolysis, for output
@@ -679,5 +546,4 @@
                iter_3d(ig,l) = iter(l)
             end do
-
 !        condensation of h2o2
 
@@ -685,16 +551,15 @@
                          ig,ptimestep,pplev,pplay,                 &
                          ztemp,zycol,dqcloud,dqscloud)
-         end if
-
-!        chemistry in upper atmosphere
-        
-         if (thermochem) then
-            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
+
+            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
 
@@ -737,5 +602,4 @@
       end do ! of do ig=1,ngrid
 
-!      stop
 !=======================================================================
 !     write outputs
@@ -753,8 +617,10 @@
             call writediagfi(ngrid,'iter','iterations',  &
                              ' ',3,iter_3d(1,1))
-            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)
+            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
             
             if (callstats) then
Index: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 2157)
+++ trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 2158)
@@ -10,5 +10,5 @@
      &                      igcm_nplus, igcm_noplus, igcm_n2plus,
      &                      igcm_hplus, igcm_hco2plus, mmol,
-     &                      igcm_he
+     &                      igcm_he, igcm_elec
       use conc_mod, only: mmean, Akknew, rnew, cpnew
       implicit none
@@ -239,4 +239,10 @@
             nbq = nbq + 1
             niq(nbq) = igcm_hco2plus
+            aki(nbq) = 0.0
+            cpi(nbq) = 0.0
+         endif
+         if(igcm_elec /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_elec
             aki(nbq) = 0.0
             cpi(nbq) = 0.0
Index: trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90	(revision 2157)
+++ trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90	(revision 2158)
@@ -106,3 +106,64 @@
          END SUBROUTINE allocate_param_iono
 
+
+!***********************************************************************
+      function temp_elect(zkm,tt,origin)
+
+!     Computes the electronic temperature, either from Viking (origin=1) 
+!     or MAVEN (origin=2) 
+
+!***********************************************************************
+      
+!     Arguments 	
+
+      real 	      tt        ! Temperature 
+      real            zkm       !  Altitude in km
+      integer         origin    ! Viking (origin=1) or MAVEN (origin=2)
+
+! local variables: 
+      real          temp_elect     ! electronic temperatures
+      real          zhanson(9),tehanson(9)
+      real          incremento
+      integer       ii, i1, i2
+
+      zhanson(1:9) = (/ 120.,130.,150.,175.,200.,225.,250.,275.,300. /)
+      tehanson(2:9) = (/ 200.,300.,500.,1250.,2000.,2200.,2400.,2500. /)
+      tehanson(1) = tt
+
+      if(origin.eq.1) then
+         if ( zkm .le. 120. ) then
+            temp_elect = tt 
+         else if(zkm .ge.300.) then
+            temp_elect=tehanson(9)
+         else
+            do ii=9,2,-1 
+               if ( zkm .lt. zhanson(ii) ) then 
+                  i1 = ii - 1
+                  i2 = ii 
+               endif
+            enddo
+            incremento=(tehanson(i2)-tehanson(i1))/(zhanson(i2)-zhanson(i1))
+            temp_elect = tehanson(i1) + (zkm-zhanson(i1)) * incremento
+	endif
+      else if(origin.eq.2) then
+         !MAVEN measured electronic temperature (Ergun et al., GRL 2015)
+         !Note that the Langmuir probe is not sensitive below ~500K, so 
+         !electronic temperatures in the lower thermosphere (<150 km) may
+         !be overestimated by this formula
+         if(zkm.le.120) then
+            temp_elect = tt
+         else 
+            temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
+         endif
+      else
+         write(*,*)'Error in function temp_elect:'
+         write(*,*)'Origin must be either 1 or 2'
+         write(*,*)'Using neutral temperature instead'
+         temp_elect = tt
+      endif
+
+      return
+
+      end function temp_elect
+
 END MODULE iono_h
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2157)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2158)
@@ -1,3 +1,3 @@
-!*****************************************************************
+!****************************************************************
 !
 !     Photochemical routine 
@@ -15,5 +15,6 @@
 subroutine photochemistry(nlayer, nq,                                   &
                           ig, lswitch, zycol, sza, ptimestep, press,    &
-                          alt, temp, dens, zmmean, dist_sol, zday,      &
+                          alt, temp, temp_elect,dens, zmmean,           &
+                          dist_sol, zday,                               &
                           surfdust1d, surfice1d, jo3, jh2o,tau, iter)
 
@@ -23,5 +24,5 @@
                            jonline
 
-use param_v4_h, only: jdistot, jdistot_b
+use param_v4_h, only: jion
 
 implicit none
@@ -42,4 +43,5 @@
 real :: alt(nlayer)           ! altitude (km)
 real :: temp(nlayer)          ! temperature (k)
+real :: temp_elect(nlayer)    ! electronic temperature (K)
 real :: dens(nlayer)          ! density (cm-3)
 real :: zmmean(nlayer)        ! mean molar mass (g/mole)
@@ -68,5 +70,5 @@
 !===================================================================
 
-integer, parameter :: nesp = 17  ! number of species in the chemical code
+integer, parameter :: nesp = 28  ! number of species in the chemical code
 
 integer :: phychemrat            ! (physical timestep)/(nominal chemical timestep)
@@ -74,5 +76,5 @@
 integer :: lswitch
 logical, save :: firstcall = .true.
-logical :: jparam                ! switch for J parameterization
+logical :: jionos                ! switch for J parameterization
 
 ! tracer indexes in the chemistry:
@@ -95,4 +97,15 @@
 integer,parameter :: i_no2  = 16
 integer,parameter :: i_n2   = 17
+integer,parameter :: i_co2plus = 18
+integer,parameter :: i_oplus   = 19
+integer,parameter :: i_o2plus  = 20
+integer,parameter :: i_noplus  = 21
+integer,parameter :: i_coplus  = 22
+integer,parameter :: i_cplus   = 23
+integer,parameter :: i_n2plus  = 24
+integer,parameter :: i_nplus   = 25
+integer,parameter :: i_hplus   = 26
+integer,parameter :: i_hco2plus= 27
+integer,parameter :: i_elec    = 28
 
 
@@ -134,5 +147,7 @@
    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_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.
 end if
@@ -145,5 +160,8 @@
                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, dens, rm, c) 
+               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, rm, c) 
 
 !===================================================================
@@ -151,5 +169,5 @@
 !===================================================================
 
-jparam= .false.
+jionos= .true.
 
 if (jonline) then
@@ -160,25 +178,70 @@
                              i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,   &
                              tau, sza, dist_sol, v_phot)
+
+      if(jionos) then
+         call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
+         do ilay=1,lswitch-1
+            call phdisrate(ig,nlayer,2,sza,ilay)
+         enddo
+         v_phot(:,14)=jion(1,:,1)
+         v_phot(:,15)=jion(1,:,2)
+         v_phot(:,16)=jion(1,:,2)
+         v_phot(:,17)=jion(1,:,3)
+         v_phot(:,18)=jion(1,:,3)
+         v_phot(:,19)=jion(1,:,4)
+         v_phot(:,20)=jion(1,:,4)
+         v_phot(:,21)=jion(2,:,1)
+         v_phot(:,22)=jion(3,:,1)
+         v_phot(:,23)=jion(10,:,1)
+         v_phot(:,24)=jion(11,:,1)
+         v_phot(:,25)=jion(11,:,2)
+         v_phot(:,26)=jion(11,:,2)
+         v_phot(:,27)=jion(8,:,1)
+         v_phot(:,28)=jion(8,:,2)
+         v_phot(:,29)=jion(8,:,2)
+         v_phot(:,30)=jion(9,:,1)
+         v_phot(:,31)=jion(12,:,1)
+      endif
+!      write(*,*)'photochemistry/205',c(:,i_co2),ig
+!      write(*,*)'photochemistry/206',v_phot(:,3),ig
+      
    else ! night
       v_phot(:,:) = 0.
    end if
-else if(jparam) then
-   call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
-   
-   do ilay=1,lswitch-1
-      call phdisrate(ig,nlayer,2,sza,ilay)
-   enddo
-   v_phot(:,1)=jdistot(2,:)
-   v_phot(:,2)=jdistot_b(2,:)
-   v_phot(:,3)=jdistot(1,:)
-   v_phot(:,4)=jdistot_b(1,:)
-   v_phot(:,5)=jdistot(7,:)
-   v_phot(:,6)=jdistot_b(7,:)
-   v_phot(:,7)=jdistot(4,:)
-   v_phot(:,8)=jdistot(6,:)
-   v_phot(:,10)=jdistot(5,:)
-   v_phot(:,11)=jdistot(10,:)
-   v_phot(:,12)=jdistot(13,:)
-   v_phot(:,13)=jdistot(8,:)
+!else if(jparam) then
+!   call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
+!   do ilay=1,lswitch-1
+!      call phdisrate(ig,nlayer,2,sza,ilay)
+!   enddo
+!   v_phot(:,1)=jdistot(2,:)
+!   v_phot(:,2)=jdistot_b(2,:)
+!   v_phot(:,3)=jdistot(1,:)
+!   v_phot(:,4)=jdistot_b(1,:)
+!   v_phot(:,5)=jdistot(7,:)
+!   v_phot(:,6)=jdistot_b(7,:)
+!   v_phot(:,7)=jdistot(4,:)
+!   v_phot(:,8)=jdistot(6,:)
+!   v_phot(:,10)=jdistot(5,:)
+!   v_phot(:,11)=jdistot(10,:)
+!   v_phot(:,12)=jdistot(13,:)
+!   v_phot(:,13)=jdistot(8,:)
+!   v_phot(:,14)=jion(1,:,1)
+!   v_phot(:,15)=jion(1,:,2)
+!   v_phot(:,16)=jion(1,:,2)
+!   v_phot(:,17)=jion(1,:,3)
+!   v_phot(:,18)=jion(1,:,3)
+!   v_phot(:,19)=jion(1,:,4)
+!   v_phot(:,20)=jion(1,:,4)
+!   v_phot(:,21)=jion(2,:,1)
+!   v_phot(:,22)=jion(3,:,1)
+!   v_phot(:,23)=jion(10,:,1)
+!   v_phot(:,24)=jion(11,:,1)
+!   v_phot(:,25)=jion(11,:,2)
+!   v_phot(:,26)=jion(11,:,2)
+!   v_phot(:,27)=jion(8,:,1)
+!   v_phot(:,28)=jion(8,:,2)
+!   v_phot(:,29)=jion(8,:,2)
+!   v_phot(:,30)=jion(9,:,1)
+!   v_phot(:,31)=jion(12,:,1)
 else
    tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
@@ -186,5 +249,4 @@
                    rm(:,i_co2), rm(:,i_o3), v_phot)
 end if
-
 ! save o3 and h2o photolysis for output
 
@@ -205,6 +267,7 @@
 hetero_ice  = .true.
 
-call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2),              &
-                   c(:,i_o), c(:,i_n2), press, temp, hetero_dust, hetero_ice, &
+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,                           &
                    surfdust1d, surfice1d, v_phot, v_3, v_4)
 
@@ -265,5 +328,4 @@
 
 !  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
-     
    cnew(:) = c(ilev,:)
 
@@ -287,4 +349,15 @@
    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
    cnew(:)   = 0.
 
@@ -305,6 +378,8 @@
                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, dens, c) 
-
+               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
 
@@ -431,5 +506,6 @@
 
  subroutine reactionrates(nlayer,                                  &
-                          lswitch, dens, co2, o2, o, n2, press, t, &
+                          lswitch, dens, co2, o2, o, n2, press, &
+                          t, t_elect,                              &
                           hetero_dust, hetero_ice,                 &
                           surfdust1d, surfice1d,                   &
@@ -448,5 +524,5 @@
 
 use comcstfi_h
-use photolysis_mod, only : nphot, nb_phot_max, &
+use photolysis_mod, only : nphot, nphotion, nb_phot_max, &
                            nb_reaction_3_max,  &
                            nb_reaction_4_max
@@ -464,4 +540,5 @@
 real, dimension(nlayer) :: press             ! pressure (hPa)
 real, dimension(nlayer) :: t                 ! temperature (K)
+real, dimension(nlayer) :: t_elect           ! electronic temperature (K)
 real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
 real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
@@ -498,4 +575,9 @@
                            d008, d009, d010, d011, d012,               &
                            e001, e002,                                 &
+                           i001, i002, i003, i004, i005, i006,         &
+                           i007, i008, i009, i010, i011, i012,         &
+                           i013, i014, i015, i016, i017, i018, i019,   &
+                           i020, i021, i022, i023, i024, i025, i026,   &
+                           i027, i028, i029,                           &
                            h001, h002, h003, h004, h005
 
@@ -504,5 +586,5 @@
 !----------------------------------------------------------------------
 
-      nb_phot       = nphot ! initialised to the number of photolysis rates
+      nb_phot       = nphot + nphotion ! initialised to the number of photolysis and photoionization rates
       nb_reaction_3 = 0
       nb_reaction_4 = 0
@@ -986,4 +1068,267 @@
       nb_reaction_4 = nb_reaction_4 + 1
       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
+      
+      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
+      
+      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(:)
 
 !----------------------------------------------------------------------
@@ -1174,4 +1519,5 @@
   mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
 
+
   loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
   loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
@@ -1187,5 +1533,7 @@
  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_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)
 
 !================================================================
@@ -1212,5 +1560,7 @@
 integer :: 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_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
 
 ! local
@@ -1334,4 +1684,146 @@
 
 !===========================================================
+!      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)
+
+!===========================================================
 !      a001 : O + O2 + CO2 -> O3 + CO2
 !===========================================================
@@ -1660,4 +2152,237 @@
 
 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)
+
+!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)
+
+!===========================================================
+!      i029 : HCO2+ + e -> OH + 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)
 
 !===========================================================
@@ -1734,5 +2459,8 @@
                            i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
                            i_n, i_n2d, i_no, i_no2, i_n2,            &
-                           dens, rm, c) 
+                           i_co2plus, i_oplus, i_o2plus, i_noplus,   &
+                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
+                           i_hplus, i_hco2plus, i_elec, dens, rm, c) 
+        
 
 !*****************************************************************
@@ -1741,5 +2469,9 @@
      &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,  &
      &                      igcm_ho2, igcm_h2o2, igcm_h2o_vap,           &
-     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
+     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
+     &                      igcm_co2plus, igcm_oplus, igcm_o2plus,       &
+     &                      igcm_noplus, igcm_coplus, igcm_cplus,        &
+     &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
+     &                      igcm_hco2plus, igcm_elec
 
       implicit none
@@ -1758,5 +2490,7 @@
       integer :: 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_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
 
       real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
@@ -1852,4 +2586,48 @@
             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
          firstcall = .false.
       end if ! of if (firstcall)
@@ -1877,4 +2655,15 @@
          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 
 
@@ -1888,5 +2677,5 @@
 
       do iesp = 1,nesp
-         do l = 1,nlayer!-1!lswitch-1
+         do l = 1,nlayer
             c(l,iesp) = rm(l,iesp)*dens(l)
          end do
@@ -1900,5 +2689,9 @@
                            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, dens, c) 
+                           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) 
+        
  
 !*****************************************************************
@@ -1907,5 +2700,9 @@
                             igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,   &
                             igcm_ho2, igcm_h2o2, igcm_h2o_vap,            &
-                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
+                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
+                            igcm_co2plus, igcm_oplus, igcm_o2plus,        &
+                            igcm_noplus, igcm_coplus, igcm_cplus,         &
+                            igcm_n2plus, igcm_nplus, igcm_hplus,          &
+                            igcm_hco2plus, igcm_elec
 
       implicit none
@@ -1923,5 +2720,8 @@
       integer :: 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_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
 
       real :: dens(nlayer)     ! total number density (molecule.cm-3) 
@@ -1962,4 +2762,15 @@
          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 
 
Index: trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90	(revision 2157)
+++ trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90	(revision 2158)
@@ -4,14 +4,14 @@
 
 ! 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 + 9 ! photolysis + quenching/heterogeneous
+  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 = 31        ! bimolecular
+  integer, parameter :: nb_reaction_4_max = 60        ! bimolecular
 
 ! spectral grid
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2157)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2158)
@@ -31,5 +31,5 @@
      &                      nuice_ref, rho_ice, rho_dust, ref_r0,
      &                      igcm_he, igcm_stormdust_mass,
-     &                      igcm_stormdust_number 
+     &                      igcm_stormdust_number
       use comsoil_h, only: inertiedat, ! soil thermal inertia
      &                     tsoil, nsoilmx ! number of subsurface layers
@@ -576,5 +576,5 @@
          if(callnlte.and.nltemodel.eq.2) call nlte_setup
          if(callnirco2.and.nircorr.eq.1) call NIR_leedat
-         if(thermochem) call chemthermos_readini
+
 
         IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
@@ -1614,5 +1614,5 @@
 c        photochemistry :
 c        --------------
-         IF (photochem .or. thermochem) then
+         IF (photochem) then
 
 !           dust and ice surface area
@@ -1665,5 +1665,5 @@
            endif
 
-         END IF  ! of IF (photochem.or.thermochem)
+         END IF  ! of IF (photochem)
 #endif
 
@@ -2396,5 +2396,5 @@
            endif ! (dustbin.ne.0)
 
-           if (thermochem .or. photochem) then
+           if (photochem) then
               do iq=1,nq
                  if (noms(iq) .ne. "dust_mass" .and.
@@ -2428,8 +2428,8 @@
      &                           (1000*mmol(iq))
 
-!                   call wstats(ngrid,"rho_"//trim(noms(iq)),
-!    $                     "Number density","cm-3",3,rhopart)
-!                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
-!    $                     "Number density","cm-3",3,rhopart)
+                   call wstats(ngrid,"rho_"//trim(noms(iq)),
+     $                   "Number density","cm-3",3,rhopart)
+                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
+     $                  "Number density","cm-3",3,rhopart)
 
 !                   vertical column (molecule.cm-2)
@@ -2461,5 +2461,5 @@
                  end if ! of if (noms(iq) .ne. "dust_mass" ...)
               end do ! of do iq=1,nq
-           end if ! of if (thermochem .or. photochem)
+           end if ! of if (photochem)
 
            end if ! of if (tracer)
@@ -2939,4 +2939,13 @@
             call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
      $           3,zdtnirco2)
+
+!            call wstats(ngrid,"q15um","15 um cooling","K/s",
+!     $           3,zdtnlte)
+!            call wstats(ngrid,"quv","UV heating","K/s",
+!     $           3,zdteuv)
+!            call wstats(ngrid,"cond","Thermal conduction","K/s",
+!     $           3,zdtconduc)
+!            call wstats(ngrid,"qnir","NIR heating","K/s",
+!     $           3,zdtnirco2)
 
          endif  !(callthermos)
