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)
