Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2006)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2007)
@@ -162,4 +162,5 @@
       real :: zlocal(nlayer)       ! Altitude (km)
       real :: zycol(nlayer,nq)     ! Composition (mole fractions)
+      real :: zmmean(nlayer)       ! Mean molecular mass (g.mole-1)
       real :: szacol               ! Solar zenith angle
       real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
@@ -169,4 +170,7 @@
       real :: em_o2(nlayer)        !  O2 nightglow emission rate      
 
+      integer :: iter(nlayer)      !  Number of chemical iterations
+                                   !  within one physical timestep 
+
 !     for output:
 
@@ -176,4 +180,6 @@
       real :: emission_no(ngrid,nlayer) !NO emission rate
       real :: emission_o2(ngrid,nlayer) !O2 emission rate
+      real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
+                                    !  within one physical timestep
 
 !=======================================================================
@@ -296,4 +302,44 @@
             niq(nbq) = i_n2
          end if
+         i_n = igcm_n
+         if (i_n == 0) then
+            if (photochem) then
+               write(*,*) "calchim: Error; no N tracer !!!"
+               stop
+            end if
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_n
+         end if
+         i_n2d = igcm_n2d
+         if (i_n2d == 0) then
+            if (photochem) then
+               write(*,*) "calchim: Error; no N2D tracer !!!"
+               stop
+            end if
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_n2d
+         end if
+         i_no = igcm_no
+         if (i_no == 0) then
+            if (photochem) then
+               write(*,*) "calchim: Error; no NO tracer !!!"
+               stop
+            end if
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_no
+         end if
+         i_no2 = igcm_no2
+         if (i_no2 == 0) then
+            if (photochem) then
+               write(*,*) "calchim: Error; no NO2 tracer !!!"
+               stop
+            end if
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_no2
+         end if
          i_h2o = igcm_h2o_vap
          if (i_h2o == 0) then
@@ -310,11 +356,8 @@
             !NO is used to determine if N chemistry is wanted
             !chemthermod=2 -> N chemistry
-            i_no = igcm_no
             if (i_no == 0) then
                write(*,*) "calchim: no NO tracer"
                write(*,*) "C/O/H themosp chemistry only "
             else
-               nbq = nbq + 1
-               niq(nbq) = i_no
                chemthermod=2
                write(*,*) "calchim: NO in traceur.def"
@@ -322,5 +365,4 @@
             end if
             ! N
-            i_n = igcm_n
             if(chemthermod == 2) then
                if (i_n == 0) then
@@ -328,50 +370,17 @@
                   write(*,*) "N is needed if NO is in traceur.def"
                   stop
-               else
-                  nbq = nbq + 1
-                  niq(nbq) = i_n
-               end if
-            else 
-               if (i_n /= 0) then
-                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
-                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
-                  stop
-               endif
-            endif    !Of if(chemthermod == 2) 
+               end if
             ! NO2
-            i_no2 = igcm_no2
-            if(chemthermod == 2) then
                if (i_no2 == 0) then
                   write(*,*) "calchim: Error; no NO2 tracer !!!"
                   write(*,*) "NO2 is needed if NO is in traceur.def"
                   stop
-               else
-                  nbq = nbq + 1
-                  niq(nbq) = i_no2
-               end if
-            else
-               if (i_no2 /= 0) then
-                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
-                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
-                  stop
-               endif
-            endif     !Of if(chemthermod == 2)
+               end if
             ! N(2D)
-            if(chemthermod == 2) then
-               i_n2d = igcm_n2d
                if (i_n2d == 0) then
                   write(*,*) "calchim: Error; no N2D !!!"
                   write(*,*) "N2D is needed if NO is in traceur.def"
                   stop
-               else
-                  nbq = nbq + 1
-                  niq(nbq) = i_n2d
-               end if
-            else
-               if (i_n2d /= 0) then
-                  write(*,*) "calchim: Error: N2D is present, but NO is not!!!"
-                  write(*,*) "Both must be in traceur.def if N chemistry wanted"
-                  stop
-               endif
+               end if
             endif    !Of if(chemthermod == 2)
             ! Ions
@@ -614,4 +623,5 @@
             zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
             zlocal(l) = zzlay(ig,l)/1000.
+            zmmean(l) = mmean(ig,l)
 
 !           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
@@ -647,8 +657,9 @@
 
          if (photochem) then
-            call photochemistry(nlayer,nq,                         &
-                                lswitch,zycol,szacol,ptimestep,    &
-                                zpress,ztemp,zdens,dist_sol,       &
-                                surfdust1d,surfice1d,jo3,taucol)
+
+            call photochemistry(nlayer,nq,                            &
+                                ig,lswitch,zycol,szacol,ptimestep,    &
+                                zpress,ztemp,zdens,zmmean,dist_sol,   &
+                                surfdust1d,surfice1d,jo3,taucol,iter)
 
 !        ozone photolysis, for output
@@ -656,4 +667,5 @@
             do l = 1,nlayer
                jo3_3d(ig,l) = jo3(l)
+               iter_3d(ig,l) = iter(l)
             end do
 
@@ -684,4 +696,5 @@
                             zu, zv, zt, zycol, ptimestep, co2ice)
          end if
+
 !=======================================================================
 !     tendencies
@@ -693,4 +706,5 @@
 
 !     tendency for the most abundant species = - sum of others
+
          do l = 1,nlayer
             iloc=maxloc(zycol(l,:))
@@ -724,4 +738,6 @@
             call writediagfi(ngrid,'jo3','j o3->o1d',    &
                              's-1',3,jo3_3d(1,1))
+            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)
@@ -729,12 +745,14 @@
                  'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
 
-           if (callstats) then
-              call wstats(ngrid,'jo3','j o3->o1d',       &
-                          's-1',3,jo3_3d(1,1))
-              call wstats(ngrid,'emission_no',           &
+            if (callstats) then
+               call wstats(ngrid,'jo3','j o3->o1d',       &
+                           's-1',3,jo3_3d(1,1))
+               call wstats(ngrid,'emission_no',           &
                    'NO nightglow emission rate','cm-3 s-1',3,emission_no)
-              call wstats(ngrid,'emission_o2',           &
+               call wstats(ngrid,'emission_o2',           &
                    'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
-           endif
+               call wstats(ngrid,'mmean','mean molecular mass',       &
+                           'g.mole-1',3,mmean(1,1))
+            endif
          end if ! of if (ngrid.gt.1)
       end if ! of if (output)
@@ -742,3 +760,2 @@
       return
       end
-
Index: unk/LMDZ.MARS/libf/aeronomars/calchim_asis.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim_asis.F90	(revision 2006)
+++ 	(revision )
@@ -1,761 +1,0 @@
-      subroutine calchim_asis(ngrid,nlayer,nq,                              &
-                         ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
-                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
-                         dqscloud,tauref,co2ice,                            &
-                         pu,pdu,pv,pdv,surfdust,surfice)
-
-      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
-                            igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2,  &
-                            igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap,   &
-                            igcm_no, igcm_n, igcm_no2, igcm_n2d,          &
-                            igcm_o2plus, igcm_co2plus, igcm_oplus,        &
-                            igcm_coplus, igcm_cplus, igcm_nplus,          &
-                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
-                            igcm_hco2plus, igcm_elec, mmol
-
-      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
-      USE comcstfi_h
-
-      implicit none
-
-!=======================================================================
-!
-!   subject:
-!   --------
-!
-!  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:
-!   ----------
-!
-!  Input:
-!
-!    ptimestep                  timestep (s)
-!    pplay(ngrid,nlayer)    Pressure at the middle of the layers (Pa)
-!    pplev(ngrid,nlayer+1)  Intermediate pressure levels (Pa)
-!    pt(ngrid,nlayer)       Temperature (K)
-!    pdt(ngrid,nlayer)      Temperature tendency (K)
-!    pu(ngrid,nlayer)       u component of the wind (ms-1)
-!    pdu(ngrid,nlayer)      u component tendency (K)
-!    pv(ngrid,nlayer)       v component of the wind (ms-1)
-!    pdv(ngrid,nlayer)      v component tendency (K)
-!    dist_sol                   distance of the sun (AU)
-!    mu0(ngrid)               cos of solar zenith angle (=1 when sun at zenith)
-!    pq(ngrid,nlayer,nq)  Advected fields, ie chemical species here
-!    pdq(ngrid,nlayer,nq) Previous tendencies on pq
-!    tauref(ngrid)            Optical depth at 7 hPa
-!    co2ice(ngrid)            co2 ice surface layer (kg.m-2)
-!    surfdust(ngrid,nlayer) dust surface area (m2/m3)
-!    surfice(ngrid,nlayer)  ice surface area (m2/m3)
-!
-!  Output:
-!
-!    dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
-!    dqschim(ngrid,nq)         ! tendencies on qsurf 
-!
-!=======================================================================
-
-#include "chimiedata.h"
-#include "callkeys.h"
-
-!     input:
-
-      integer,intent(in) :: ngrid    ! number of atmospheric columns
-      integer,intent(in) :: nlayer   ! number of atmospheric layers
-      integer,intent(in) :: nq       ! number of tracers
-      real :: ptimestep
-      real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
-      real :: zzlay(ngrid,nlayer)    ! pressure at the middle of the layers
-      real :: pplev(ngrid,nlayer+1)  ! intermediate pressure levels
-      real :: zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries
-      real :: pt(ngrid,nlayer)       ! temperature
-      real :: pdt(ngrid,nlayer)      ! temperature tendency
-      real :: pu(ngrid,nlayer)       ! u component of the wind (m.s-1)
-      real :: pdu(ngrid,nlayer)      ! u component tendency
-      real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
-      real :: pdv(ngrid,nlayer)      ! v component tendency
-      real :: dist_sol               ! distance of the sun (AU)
-      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
-      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
-      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
-      real :: zday                   ! date (time since Ls=0, in martian days)
-      real :: tauref(ngrid)          ! optical depth at 7 hPa
-      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
-      real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
-      real :: surfice(ngrid,nlayer)  !  ice surface area (m2/m3)
-
-!     output:
-
-      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
-      real :: dqschim(ngrid,nq)         ! tendencies on qsurf 
-      real :: dqcloud(ngrid,nlayer,nq)  ! tendencies on pq due to condensation
-      real :: dqscloud(ngrid,nq)        ! tendencies on qsurf 
-
-!     local variables:
-
-      integer,save :: nbq                   ! number of tracers used in the chemistry
-      integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
-      integer :: iloc(1)                    ! index of major species
-      integer :: ig,l,i,iq,iqmax
-      integer :: foundswitch, lswitch
-      integer,save :: chemthermod
-
-      integer,save :: i_co2  = 0
-      integer,save :: i_co   = 0
-      integer,save :: i_o    = 0
-      integer,save :: i_o1d  = 0
-      integer,save :: i_o2   = 0
-      integer,save :: i_o3   = 0
-      integer,save :: i_h    = 0
-      integer,save :: i_h2   = 0
-      integer,save :: i_oh   = 0
-      integer,save :: i_ho2  = 0
-      integer,save :: i_h2o2 = 0
-      integer,save :: i_ch4  = 0
-      integer,save :: i_n2   = 0
-      integer,save :: i_h2o  = 0
-      integer,save :: i_n    = 0
-      integer,save :: i_no   = 0
-      integer,save :: i_no2  = 0
-      integer,save :: i_n2d  = 0
-      integer,save :: i_co2plus=0
-      integer,save :: i_oplus=0
-      integer,save :: i_o2plus=0
-      integer,save :: i_coplus=0
-      integer,save :: i_cplus=0
-      integer,save :: i_nplus=0
-      integer,save :: i_noplus=0
-      integer,save :: i_n2plus=0
-      integer,save :: i_hplus=0
-      integer,save :: i_hco2plus=0
-      integer,save :: i_elec=0
-
-      integer :: ig_vl1
-
-      real    :: latvl1, lonvl1
-      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
-                                       ! new mole fraction after
-      real    :: zt(ngrid,nlayer)      ! temperature
-      real    :: zu(ngrid,nlayer)      ! u component of the wind
-      real    :: zv(ngrid,nlayer)      ! v component of the wind
-      real    :: taucol                ! optical depth at 7 hPa
-
-      logical,save :: firstcall = .true.
-      logical,save :: depos = .false.  ! switch for dry deposition
-
-!     for each column of atmosphere:
-
-      real :: zpress(nlayer)       ! Pressure (mbar)
-      real :: zdens(nlayer)        ! Density  (cm-3)
-      real :: ztemp(nlayer)        ! Temperature (K)
-      real :: zlocal(nlayer)       ! Altitude (km)
-      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
-      real :: zmmean(nlayer)       ! Mean molecular mass (g.mole-1)
-      real :: szacol               ! Solar zenith angle
-      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
-      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
-      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (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 
-
-!     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 :: emission_no(ngrid,nlayer) !NO emission rate
-      real :: emission_o2(ngrid,nlayer) !O2 emission rate
-      real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
-                                    !  within one physical timestep
-
-!=======================================================================
-!     initialization of the chemistry (first call only)
-!=======================================================================
-
-      if (firstcall) then
-
-         if (photochem) then
-            print*,'calchim: Read photolysis lookup table'
-            call read_phototable
-         end if
-         ! 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.
-         i_co2 = igcm_co2
-         if (i_co2 == 0) then
-            write(*,*) "calchim: Error; no CO2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_co2
-         end if
-         i_co = igcm_co
-         if (i_co == 0) then
-            write(*,*) "calchim: Error; no CO tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_co
-         end if
-         i_o = igcm_o
-         if (i_o == 0) then
-            write(*,*) "calchim: Error; no O tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o
-         end if
-         i_o1d = igcm_o1d
-         if (i_o1d == 0) then
-            write(*,*) "calchim: Error; no O1D tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o1d
-         end if
-         i_o2 = igcm_o2
-         if (i_o2 == 0) then
-            write(*,*) "calchim: Error; no O2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o2
-         end if
-         i_o3 = igcm_o3
-         if (i_o3 == 0) then
-            write(*,*) "calchim: Error; no O3 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o3
-         end if
-         i_h = igcm_h
-         if (i_h == 0) then
-            write(*,*) "calchim: Error; no H tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h
-         end if
-         i_h2 = igcm_h2
-         if (i_h2 == 0) then
-            write(*,*) "calchim: Error; no H2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h2
-         end if
-         i_oh = igcm_oh
-         if (i_oh == 0) then
-            write(*,*) "calchim: Error; no OH tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_oh
-         end if
-         i_ho2 = igcm_ho2
-         if (i_ho2 == 0) then
-            write(*,*) "calchim: Error; no HO2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_ho2
-         end if
-         i_h2o2 = igcm_h2o2
-         if (i_h2o2 == 0) then
-            write(*,*) "calchim: Error; no H2O2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h2o2
-         end if
-         i_ch4 = igcm_ch4
-         if (i_ch4 == 0) then
-            write(*,*) "calchim: Error; no CH4 tracer !!!"
-            write(*,*) "CH4 will be ignored in the chemistry"
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_ch4
-         end if
-         i_n2 = igcm_n2
-         if (i_n2 == 0) then
-            write(*,*) "calchim: Error; no N2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n2
-         end if
-         i_n = igcm_n
-         if (i_n == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no N tracer !!!"
-               stop
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n
-         end if
-         i_n2d = igcm_n2d
-         if (i_n2d == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no N2D tracer !!!"
-               stop
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n2d
-         end if
-         i_no = igcm_no
-         if (i_no == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no NO tracer !!!"
-               stop
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_no
-         end if
-         i_no2 = igcm_no2
-         if (i_no2 == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no NO2 tracer !!!"
-               stop
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_no2
-         end if
-         i_h2o = igcm_h2o_vap
-         if (i_h2o == 0) then
-            write(*,*) "calchim: Error; no water vapor tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            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
-
-         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
-               
-         firstcall = .false.
-      end if ! if (firstcall)
-
-! Initializations
-
-      zycol(:,:)    = 0.
-      dqchim(:,:,:) = 0.
-      dqschim(:,:)  = 0.
-
-!     latvl1= 22.27
-!     lonvl1= -47.94
-!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +    &
-!             int(1.5+(lonvl1+180)*iim/360.)
-
-!=======================================================================
-!     loop over grid
-!=======================================================================
-
-      do ig = 1,ngrid
-         
-         foundswitch = 0
-         do l = 1,nlayer
-            do i = 1,nbq
-               iq = niq(i) ! get tracer index
-               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-            end do
-            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
-            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
-            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
-            zpress(l) = pplay(ig,l)/100.
-            ztemp(l)  = zt(ig,l)
-            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
-            zlocal(l) = zzlay(ig,l)/1000.
-            zmmean(l) = mmean(ig,l)
-
-!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
-
-            surfdust1d(l) = surfdust(ig,l)*1.e-2
-            surfice1d(l)  = surfice(ig,l)*1.e-2
-
-!           search for switch index between regions
-
-            if (photochem .and. thermochem) then
-               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
-
-         end do ! of do l=1,nlayer
-
-         szacol = acos(mu0(ig))*180./pi
-         taucol = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
-
-!=======================================================================
-!     call chemical subroutines
-!=======================================================================
-
-!        chemistry in lower atmosphere
-
-         if (photochem) then
-
-            call photochemistry_asis(nlayer,nq,                       &
-                                ig,lswitch,zycol,szacol,ptimestep,    &
-                                zpress,ztemp,zdens,zmmean,dist_sol,   &
-                                surfdust1d,surfice1d,jo3,taucol,iter)
-
-!        ozone photolysis, for output
-
-            do l = 1,nlayer
-               jo3_3d(ig,l) = jo3(l)
-               iter_3d(ig,l) = iter(l)
-            end do
-
-!        condensation of h2o2
-
-            call perosat(ngrid, nlayer, nq,                        &
-                         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
-         end if
-
-!        dry deposition
-
-         if (depos) then
-            call deposition(ngrid, nlayer, nq,                      &
-                            ig, ig_vl1, pplay, pplev, zzlay, zzlev, & 
-                            zu, zv, zt, zycol, ptimestep, co2ice)
-         end if
-
-!=======================================================================
-!     tendencies
-!=======================================================================
-
-!     index of the most abundant species at each level
-
-!         major(:) = maxloc(zycol, dim = 2)
-
-!     tendency for the most abundant species = - sum of others
-
-         do l = 1,nlayer
-            iloc=maxloc(zycol(l,:))
-            iqmax=iloc(1)
-            do i = 1,nbq
-               iq = niq(i) ! get tracer index
-               if (iq /= iqmax) then
-                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)  &
-                                   - zq(ig,l,iq))/ptimestep
-                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax)              &
-                                     - dqchim(ig,l,iq) 
-               end if
-            end do
-         end do ! of do l = 1,nlayer
-
-!=======================================================================
-!     end of loop over grid
-!=======================================================================
-
-      end do ! of do ig=1,ngrid
-
-!=======================================================================
-!     write outputs
-!=======================================================================
-
-! value of parameter 'output' to trigger writting of outputs
-! is set above at the declaration of the variable.
-
-      if (photochem .and. output) then
-         if (ngrid > 1) then
-            call writediagfi(ngrid,'jo3','j o3->o1d',    &
-                             's-1',3,jo3_3d(1,1))
-            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 (callstats) then
-               call wstats(ngrid,'jo3','j o3->o1d',       &
-                           's-1',3,jo3_3d(1,1))
-               call wstats(ngrid,'emission_no',           &
-                   'NO nightglow emission rate','cm-3 s-1',3,emission_no)
-               call wstats(ngrid,'emission_o2',           &
-                   'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
-               call wstats(ngrid,'mmean','mean molecular mass',       &
-                           'g.mole-1',3,mmean(1,1))
-            endif
-         end if ! of if (ngrid.gt.1)
-      end if ! of if (output)
-
-      return
-      end
Index: unk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 2006)
+++ 	(revision )
@@ -1,1570 +1,0 @@
-c*****************************************************************
-c
-c     Photochemical routine 
-c
-c     Author: Franck Lefevre
-c     ------
-c
-c     Version: 11/12/2013
-c
-c*****************************************************************
-
-      subroutine photochemistry(nlayer, nq,
-     $                  lswitch, zycol, sza, ptimestep, 
-     $                  press,temp, dens, dist_sol, surfdust1d,
-     $                  surfice1d, jo3, tau)
-
-      implicit none
-
-#include "chimiedata.h" 
-#include "callkeys.h"
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer, intent(in) :: nq     ! number of tracers
-      
-      real :: sza                   ! solar zenith angle (deg)
-      real :: ptimestep             ! physics timestep (s)
-      real :: press(nlayer)         ! pressure (hpa)
-      real :: temp(nlayer)          ! temperature (k)
-      real :: dens(nlayer)          ! density (cm-3)
-      real :: dist_sol              ! sun distance (au) 
-      real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
-      real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
-      real :: tau                   ! optical depth at 7 hpa
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     input/output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      
-      real :: zycol(nlayer,nq)      ! chemical species volume mixing ratio
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      
-      real :: jo3(nlayer)  ! photodissociation rate o3 -> o1d
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer, parameter :: nesp = 16 ! number of species in the chemistry
-
-      integer :: phychemrat       ! ratio physics timestep/chemistry timestep
-      integer :: istep
-      integer :: i_co2,i_o3,j_o3_o1d,lev
-      integer :: lswitch
-
-      real :: stimestep           ! standard timestep for the chemistry (s)
-      real :: ctimestep           ! real timestep for the chemistry (s) 
-      real :: rm(nlayer,nesp)     ! species volume mixing ratio 
-      real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
-      real :: rmco2(nlayer)       ! co2 mixing ratio
-      real :: rmo3(nlayer)        ! ozone mixing ratio
-
-c     reaction rates
-
-      real :: a001(nlayer), a002(nlayer), a003(nlayer)
-      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
-     $        b004(nlayer), b005(nlayer), b006(nlayer),
-     $        b007(nlayer), b008(nlayer), b009(nlayer)
-      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
-     $        c004(nlayer), c005(nlayer), c006(nlayer), 
-     $        c007(nlayer), c008(nlayer), c009(nlayer), 
-     $        c010(nlayer), c011(nlayer), c012(nlayer),
-     $        c013(nlayer), c014(nlayer), c015(nlayer), 
-     $        c016(nlayer), c017(nlayer), c018(nlayer)
-      real :: d001(nlayer), d002(nlayer), d003(nlayer)
-      real :: e001(nlayer), e002(nlayer), e003(nlayer)
-      real :: h001(nlayer), h002(nlayer), h003(nlayer),
-     $        h004(nlayer), h005(nlayer)
-      real :: t001(nlayer), t002(nlayer), t003(nlayer)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     stimestep  : standard timestep for the chemistry (s)       c
-c     ctimestep  : real timestep for the chemistry (s)           c
-c     phychemrat : ratio physical/chemical timestep              c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      stimestep = 600. ! standard value : 10 mn
-
-      phychemrat = nint(ptimestep/stimestep)
-
-      ctimestep = ptimestep/real(phychemrat)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialisation of chemical species and families            c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      call gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm) 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     compute reaction rates                                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      call chemrates(nlayer,
-     $               lswitch, dens, press, temp,
-     $               surfdust1d, surfice1d,
-     $               a001, a002, a003,
-     $               b001, b002, b003, b004, b005, b006,
-     $               b007, b008, b009,
-     $               c001, c002, c003, c004, c005, c006,
-     $               c007, c008, c009, c010, c011, c012,
-     $               c013, c014, c015, c016, c017, c018,
-     $               d001, d002, d003,
-     $               e001, e002, e003,
-     $               h001, h002, h003, h004, h005,
-     $               t001, t002, t003, tau)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     interpolation of photolysis rates in the lookup table      c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      i_co2 = 1
-      i_o3  = 6 
-
-      do lev = 1,lswitch-1
-         rmco2(lev) = rm(lev,i_co2)
-         rmo3(lev)  = rm(lev, i_o3)
-      end do
-
-      call photolysis(nlayer, lswitch, press, temp, sza, tau, dist_sol,
-     $                rmco2, rmo3, j)
-
-      j_o3_o1d = 5
-
-      do lev = 1,lswitch-1
-         jo3(lev) = j(lev,j_o3_o1d)
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over time within the physical timestep                c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do istep = 1,phychemrat
-         call chimie(nlayer, lswitch, nesp, rm, j, dens, ctimestep, 
-     $               press, temp, sza, dist_sol,
-     $               a001, a002, a003,
-     $               b001, b002, b003, b004, b005, b006,
-     $               b007, b008, b009,
-     $               c001, c002, c003, c004, c005, c006,
-     $               c007, c008, c009, c010, c011, c012,
-     $               c013, c014, c015, c016, c017, c018,
-     $               d001, d002, d003,
-     $               e001, e002, e003,
-     $               h001, h002, h003, h004, h005,
-     $               t001, t002, t003)
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     save chemical species for the gcm                          c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      call chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm) 
-
-      contains
-
-
-c*****************************************************************
-
-      subroutine chimie(nlayer,
-     $                  lswitch, nesp, rm, j, dens, dt, 
-     $                  press, t, sza, dist_sol,
-     $                  a001, a002, a003,
-     $                  b001, b002, b003, b004, b005, b006,
-     $                  b007, b008, b009,
-     $                  c001, c002, c003, c004, c005, c006,
-     $                  c007, c008, c009, c010, c011, c012,
-     $                  c013, c014, c015, c016, c017, c018,
-     $                  d001, d002, d003,
-     $                  e001, e002, e003,
-     $                  h001, h002, h003, h004, h005,
-     $                  t001, t002, t003)
-
-c*****************************************************************
-
-      use tracer_mod, only: igcm_ch4
-      implicit none 
-
-#include "chimiedata.h" 
-#include "callkeys.h" 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer :: lswitch            ! interface level between chemistries
-      integer :: nesp               ! number of species
-
-      real :: dens(nlayer)          ! density (cm-3)
-      real :: dt                    ! chemistry timestep (s) 
-      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
-      real :: press(nlayer)         ! pressure (hpa)
-      real :: t(nlayer)             ! temperature (k)
-      real :: sza                   ! solar zenith angle (deg)
-      real :: dist_sol              ! sun distance (au) 
-
-c     reaction rates
-
-      real :: a001(nlayer), a002(nlayer), a003(nlayer)
-      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
-     $        b004(nlayer), b005(nlayer), b006(nlayer),
-     $        b007(nlayer), b008(nlayer), b009(nlayer)
-      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
-     $        c004(nlayer), c005(nlayer), c006(nlayer), 
-     $        c007(nlayer), c008(nlayer), c009(nlayer), 
-     $        c010(nlayer), c011(nlayer), c012(nlayer),
-     $        c013(nlayer), c014(nlayer), c015(nlayer), 
-     $        c016(nlayer), c017(nlayer), c018(nlayer)
-      real :: d001(nlayer), d002(nlayer), d003(nlayer)
-      real :: e001(nlayer), e002(nlayer), e003(nlayer)
-      real :: h001(nlayer), h002(nlayer), h003(nlayer),
-     $        h004(nlayer), h005(nlayer)
-      real :: t001(nlayer), t002(nlayer), t003(nlayer)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     input/output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: rm(nlayer,nesp)  ! volume mixing ratios 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real, parameter :: hetero_ice  = 1. ! switch for het. chem. on ice clouds
-      real, parameter :: hetero_dust = 0. ! switch for het. chem. on dust
-
-      integer :: iesp, iter, l
-      integer :: i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
-     $           i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
-      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
-     $           j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
-     $           j_ho2, j_no, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
-     $           j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
-     $           j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
-     $           j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
-     $           j_hno3, j_hno4 
-
-      integer, parameter :: niter = 5        ! iterations in the chemical scheme
-
-      real :: cc0(nlayer,nesp)        ! initial number density (cm-3)
-      real :: cc(nlayer,nesp)         ! final number density (cm-3)
-      real :: nox(nlayer)             ! nox number density (cm-3)
-      real :: no(nlayer)              ! no  number density (cm-3)
-      real :: no2(nlayer)             ! no2 number density (cm-3)
-      real :: production(nlayer,nesp) ! production rate
-      real :: loss(nlayer,nesp)       ! loss rate
-
-      real :: ro_o3, rh_ho2, roh_ho2, rno2_no
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     tracer numbering in the chemistry
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      i_co2    =  1
-      i_co     =  2
-      i_o      =  3
-      i_o1d    =  4
-      i_o2     =  5
-      i_o3     =  6
-      i_h      =  7
-      i_h2     =  8
-      i_oh     =  9
-      i_ho2    = 10 
-      i_h2o2   = 11
-      i_ch4    = 12
-      i_h2o    = 13
-      i_n2     = 14
-      i_hox    = 15
-      i_ox     = 16
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     numbering of photolysis rates
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      j_o2_o         =  1      ! o2 + hv     -> o + o
-      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
-      j_co2_o        =  3      ! co2 + hv    -> co + o
-      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
-      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
-      j_o3_o         =  6      ! o3 + hv     -> o2 + o
-      j_h2o          =  7      ! h2o + hv    -> h + oh
-      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
-      j_ho2          =  9      ! ho2 + hv    -> oh + o
-      j_no           =  10     ! no + hv     -> n + o
-      j_no2          =  11     ! no2 + hv    -> no + o
-      j_hno3         =  12     ! hno3 + hv   -> oh + no2
-      j_hno4         =  13     ! hno4 + hv   -> ho2 + no2
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     volume mixing ratio -> density conversion
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do iesp = 1,nesp
-         do l = 1,lswitch-1
-            cc0(l,iesp) = rm(l,iesp)*dens(l)
-            cc(l,iesp)  = cc0(l,iesp)
-         end do
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     co2 and nox number densities (cm-3)    
-c      
-c     nox mixing ratio: 6.e-10 (yung and demore, 1999)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         nox(l) = 6.e-10*dens(l)
-      end do 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over iterations in the chemical scheme
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do iter = 1,niter
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        nox species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no2/no
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-
-         rno2_no = (d002(l)*cc(l,i_o3) + d003(l)*cc(l,i_ho2))
-     $            /(j(l,j_no2) +
-     $              d001(l)*max(cc(l,i_o),1.e-30*dens(l)))
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         no(l) = nox(l)/(1. + rno2_no)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         no2(l) = no(l)*rno2_no
-
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        hox species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        photochemical equilibrium for oh and ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h/ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-
-         rh_ho2 = (c001(l)*cc(l,i_o)
-     $           + c004(l)*cc(l,i_h)
-     $           + c005(l)*cc(l,i_h)
-     $           + c006(l)*cc(l,i_h)
-     $           + c007(l)*cc(l,i_oh)
-     $           + 2.*c008(l)*cc(l,i_ho2)
-     $           + c015(l)*cc(l,i_o3)
-     $           + 2.*c016(l)*cc(l,i_ho2)
-     $           + d003(l)*no(l)             ! ajout 20090401
-     $           + j(l,j_ho2)
-     $           + h001(l)*hetero_ice 
-     $           + h003(l)*hetero_dust)
-     $          /( c011(l)*cc(l,i_o2)
-     $           + t001(l)*cc(l,i_h2o)
-     $            /max(cc(l,i_h),dens(l)*1.e-30)   ! ajout 20090401
-     $           )
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oh/ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         roh_ho2 = (c001(l)*cc(l,i_o)
-     $            + c003(l)*cc(l,i_o3)*rh_ho2
-     $            + 2.*c004(l)*cc(l,i_h)
-     $            + 2.*c008(l)*cc(l,i_ho2)
-     $            + c015(l)*cc(l,i_o3)
-     $            + d003(l)*no(l)
-     $            + j(l,j_ho2)
-     $            + 2.*b002(l)*cc(l,i_o1d)*cc(l,i_h2o)         ! ajout 20101210
-     $             /max(cc(l,i_ho2),dens(l)*1.e-30)            ! ajout 20101210
-     $            + b003(l)*cc(l,i_o1d)*cc(l,i_h2)             ! ajout 20101210
-     $             /max(cc(l,i_ho2),dens(l)*1.e-30)            ! ajout 20101210
-     $            + j(l,j_h2o)*cc(l,i_h2o)
-     $             /max(cc(l,i_ho2),dens(l)*1.e-30)
-     $            + t001(l)*cc(l,i_h2o)               ! suppression 20090401
-     $             /max(cc(l,i_ho2),dens(l)*1.e-30)   ! suppression 20090401
-     $            )
-     $            /(c002(l)*cc(l,i_o)
-     $            + c007(l)*cc(l,i_ho2)
-     $            + c009(l)*cc(l,i_h2o2)         ! ajout 20090401
-     $            + 2.*c013(l)*cc(l,i_oh)        ! ajout 20090401
-     $            + 2.*c017(l)*cc(l,i_oh)        ! ajout 20090401
-     $            + e001(l)*cc(l,i_co)
-     $            + h002(l)*hetero_ice)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         cc(l,i_h) = cc(l,i_hox)
-     $                 /(1. + (1. + roh_ho2)/rh_ho2)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         cc(l,i_ho2) = cc(l,i_h)/rh_ho2
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oh
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
-
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ox species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        day:
-c        - o1d at photochemical equilibrium 
-c        - o3 at photochemical equilibrium
-c        - continuity equation for ox
-c        night:
-c        - o1d = 0
-c        - continuity equation for o3
-c        - continuity equation for o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-        
-      if (sza .le. 95.) then
-
-         do l = 1,lswitch-1
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o(1d)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            cc(l,i_o1d) = (j(l,j_co2_o1d)*cc(l,i_co2)
-     $                   + j(l,j_o2_o1d)*cc(l,i_o2)
-     $                   + j(l,j_o3_o1d)*cc(l,i_o3))
-     $                   /(b001(l)*cc(l,i_co2)
-     $                   + b002(l)*cc(l,i_h2o)
-     $                   + b003(l)*cc(l,i_h2)
-     $                   + b004(l)*cc(l,i_o2)
-     $                   + b005(l)*cc(l,i_o3)
-     $                   + b006(l)*cc(l,i_o3))
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o/o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            ro_o3 = (j(l,j_o3_o1d) + j(l,j_o3_o)
-     $              + a003(l)*cc(l,i_o)
-     $              + c003(l)*cc(l,i_h)
-     $              + c014(l)*cc(l,i_oh)
-     $              + c015(l)*cc(l,i_ho2)
-     $              )
-     $              /(a001(l)*cc(l,i_o2))
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            cc(l,i_o) = cc(l,i_o3)*ro_o3
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ox = o + o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            production(l,i_ox) =
-     $                   + j(l,j_co2_o)*cc(l,i_co2)
-     $                   + j(l,j_co2_o1d)*cc(l,i_co2)
-     $                   + j(l,j_ho2)*cc(l,i_ho2)
-     $                   + 2.*j(l,j_o2_o)*cc(l,i_o2)
-     $                   + 2.*j(l,j_o2_o1d)*cc(l,i_o2)
-     $                   + c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                   + c013(l)*cc(l,i_oh)*cc(l,i_oh)
-     $                   + d003(l)*cc(l,i_ho2)*no(l)
-
-            loss(l,i_ox) = 2.*a002(l)*cc(l,i_o)*cc(l,i_o)
-     $                   + 2.*a003(l)*cc(l,i_o)*cc(l,i_o3)
-     $                   + c001(l)*cc(l,i_ho2)*cc(l,i_o)
-     $                   + c002(l)*cc(l,i_oh)*cc(l,i_o)
-     $                   + c003(l)*cc(l,i_h)*cc(l,i_o3)
-     $                   + c012(l)*cc(l,i_o)*cc(l,i_h2o2)
-     $                   + c014(l)*cc(l,i_o3)*cc(l,i_oh)
-     $                   + c015(l)*cc(l,i_o3)*cc(l,i_ho2)
-     $                   + d001(l)*cc(l,i_o)*no2(l)
-     $                   + e002(l)*cc(l,i_o)*cc(l,i_co)
-
-            loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
-
-         end do
-
-      else
-
-         do l = 1,lswitch-1
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o(1d)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            cc(l,i_o1d) = 0.
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            production(l,i_o3) = a001(l)*cc(l,i_o2)*cc(l,i_o)
-
-            loss(l,i_o3) = a003(l)*cc(l,i_o)
-     $                   + c003(l)*cc(l,i_h)
-     $                   + c014(l)*cc(l,i_oh)
-     $                   + c015(l)*cc(l,i_ho2)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-            production(l,i_o) = c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                        + c013(l)*cc(l,i_oh)*cc(l,i_oh)
-
-            loss(l,i_o)  =  a001(l)*cc(l,i_o2)
-     $                   + 2.*a002(l)*cc(l,i_o)
-     $                   + a003(l)*cc(l,i_o3)
-     $                   + c001(l)*cc(l,i_ho2)
-     $                   + c002(l)*cc(l,i_oh)
-     $                   + c012(l)*cc(l,i_h2o2)
-     $                   + e002(l)*cc(l,i_co)
-
-         end do
-      end if
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     other species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        co2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_co2) = e001(l)*cc(l,i_oh)*cc(l,i_co)
-     $                       + e002(l)*cc(l,i_o)*cc(l,i_co)
-     $                       + t002(l)*cc(l,i_ch4)*16./44. ! ajout 20090401
-     $                       + t003(l)*cc(l,i_co2)*8./44.  ! ajout 20090401
-
-         loss(l,i_co2) = j(l,j_co2_o)
-     $                 + j(l,j_co2_o1d)
-     $                 + t003(l)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        co
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_co) = j(l,j_co2_o)*cc(l,i_co2)
-     $                      + j(l,j_co2_o1d)*cc(l,i_co2)
-     $                      + t003(l)*cc(l,i_co2)
-
-         loss(l,i_co) = e001(l)*cc(l,i_oh)
-     $                + e002(l)*cc(l,i_o)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_o2) = 
-     $                  j(l,j_o3_o)*cc(l,i_o3)
-     $                + j(l,j_o3_o1d)*cc(l,i_o3)
-     $                + a002(l)*cc(l,i_o)*cc(l,i_o)
-     $                + 2.*a003(l)*cc(l,i_o)*cc(l,i_o3)
-     $                + 2.*b005(l)*cc(l,i_o1d)*cc(l,i_o3)
-     $                + b006(l)*cc(l,i_o1d)*cc(l,i_o3)
-     $                + c001(l)*cc(l,i_o)*cc(l,i_ho2)
-     $                + c002(l)*cc(l,i_o)*cc(l,i_oh)
-     $                + c003(l)*cc(l,i_h)*cc(l,i_o3)
-     $                + c005(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                + c007(l)*cc(l,i_oh)*cc(l,i_ho2)
-     $                + c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                + c014(l)*cc(l,i_o3)*cc(l,i_oh)
-     $                + 2.*c015(l)*cc(l,i_o3)*cc(l,i_ho2)
-     $                + c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                + d001(l)*cc(l,i_o)*no2(l)
-
-         loss(l,i_o2) = j(l,j_o2_o)
-     $                + j(l,j_o2_o1d)
-     $                + a001(l)*cc(l,i_o)
-     $                + c011(l)*cc(l,i_h)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_h2) = c005(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                      + c018(l)*cc(l,i_h)*cc(l,i_h)
-
-         loss(l,i_h2) = b003(l)*cc(l,i_o1d)
-     $                + c010(l)*cc(l,i_oh)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_h2o) = 
-     $                   c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                 + c007(l)*cc(l,i_oh)*cc(l,i_ho2)
-     $                 + c009(l)*cc(l,i_oh)*cc(l,i_h2o2)
-     $                 + c010(l)*cc(l,i_oh)*cc(l,i_h2)
-     $                 + c013(l)*cc(l,i_oh)*cc(l,i_oh)
-     $                 + h004(l)*cc(l,i_h2o2)*hetero_ice
-
-         loss(l,i_h2o) = j(l,j_h2o)
-     $                 + b002(l)*cc(l,i_o1d)
-     $                 + t001(l)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2o2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_h2o2) =
-     $                    c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                  + c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                  + c017(l)*cc(l,i_oh)*cc(l,i_oh)
-c    $                  + 0.5*h001(l)*cc(l,i_ho2)*hetero_ice
-c    $                  + 0.5*h002(l)*cc(l,i_oh)*hetero_ice
-
-         loss(l,i_h2o2) = j(l,j_h2o2)
-     $                  + c009(l)*cc(l,i_oh)
-     $                  + c012(l)*cc(l,i_o)
-     $                  + h004(l)*hetero_ice
-     $                  + h005(l)*hetero_dust
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        hox = h + oh + ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_hox) = 
-     $                   2.*j(l,j_h2o)*cc(l,i_h2o)
-     $                 + 2.*j(l,j_h2o2)*cc(l,i_h2o2)
-     $                 + 2.*b002(l)*cc(l,i_o1d)*cc(l,i_h2o)
-     $                 + 2.*b003(l)*cc(l,i_o1d)*cc(l,i_h2)
-     $                 + 2.*c012(l)*cc(l,i_o)*cc(l,i_h2o2)
-     $                 + 2.*t001(l)*cc(l,i_h2o)
-
-         loss(l,i_hox) = 2.*c005(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                 + 2.*c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                 + 2.*c007(l)*cc(l,i_oh)*cc(l,i_ho2)
-     $                 + 2.*c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                 + 2.*c013(l)*cc(l,i_oh)*cc(l,i_oh)
-     $                 + 2.*c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
-     $                 + 2.*c017(l)*cc(l,i_oh)*cc(l,i_oh)
-     $                 + 2.*c018(l)*cc(l,i_h)*cc(l,i_h)
-     $                 + h001(l)*cc(l,i_ho2)*hetero_ice
-     $                 + h002(l)*cc(l,i_oh)*hetero_ice
-     $                 + h003(l)*cc(l,i_ho2)*hetero_dust
-
-         loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ch4
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         production(l,i_ch4) = 0.
-
-         if (igcm_ch4/=0) then
-           loss(l,i_ch4) = j(l,j_ch4_ch3_h)
-     $                   + j(l,j_ch4_1ch2_h2)
-     $                   + j(l,j_ch4_3ch2_h_h)
-     $                   + j(l,j_ch4_ch_h2_h)
-     $                   + b007(l)*cc(l,i_o1d) 
-     $                   + b008(l)*cc(l,i_o1d)
-     $                   + b009(l)*cc(l,i_o1d)
-     $                   + e003(l)*cc(l,i_oh)
-         else
-           loss(l,i_ch4) = 0.
-         endif
-
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     update number densities
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-c     long-lived species
-
-      do l = 1,lswitch-1
-         cc(l,i_co2) = (cc0(l,i_co2) + production(l,i_co2)*dt)
-     $              /(1. + loss(l,i_co2)*dt)
-         cc(l,i_co)  = (cc0(l,i_co) + production(l,i_co)*dt)
-     $              /(1. + loss(l,i_co)*dt)
-         cc(l,i_o2)  = (cc0(l,i_o2) + production(l,i_o2)*dt)
-     $              /(1. + loss(l,i_o2)*dt)
-         cc(l,i_h2)  = (cc0(l,i_h2) + production(l,i_h2)*dt)
-     $              /(1. + loss(l,i_h2)*dt)
-         cc(l,i_h2o2)= (cc0(l,i_h2o2) + production(l,i_h2o2)*dt)
-     $              /(1. + loss(l,i_h2o2)*dt)
-         cc(l,i_h2o) = (cc0(l,i_h2o) + production(l,i_h2o)*dt)
-     $              /(1. + loss(l,i_h2o)*dt)
-         cc(l,i_hox) = (cc0(l,i_hox) + production(l,i_hox)*dt)
-     $              /(1. + loss(l,i_hox)*dt)
-         cc(l,i_ch4) = (cc0(l,i_ch4) + production(l,i_ch4)*dt)
-     $              /(1. + loss(l,i_ch4)*dt)
-      end do
-
-c     ox species
-
-      if (sza .le. 95.) then
-         do l = 1,lswitch-1
-            cc(l,i_ox) = (cc0(l,i_ox) + production(l,i_ox)*dt)
-     $                   /(1. + loss(l,i_ox)*dt)
-         end do
-      else
-         do l = 1,lswitch-1
-            cc(l,i_o)  = (cc0(l,i_o) + production(l,i_o)*dt)
-     $                   /(1. + loss(l,i_o)*dt)
-            cc(l,i_o3) = (cc0(l,i_o3) + production(l,i_o3)*dt)
-     $                   /(1. + loss(l,i_o3)*dt)
-            cc(l,i_ox) = cc(l,i_o) + cc(l,i_o3)
-         end do
-      end if 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     end of loop over iterations
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      end do
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     density -> volume mixing ratio conversion
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do iesp = 1,nesp
-         do l = 1,lswitch-1
-            rm(l,iesp) = max(cc(l,iesp)/dens(l), 1.e-30)
-         end do
-      end do
-
-      end subroutine chimie
-
-c*****************************************************************
-
-      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
-
-c*****************************************************************
-
-      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
-     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
-     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
-     &                      igcm_ch4
-
-      implicit none
-
-#include "callkeys.h"
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer, intent(in) :: nq     ! number of tracers
-      real :: zycol(nlayer,nq)      ! species volume mixing ratio in the gcm
-
-      integer :: nesp               ! number of species in the chemistry
-      integer :: lswitch            ! interface level between chemistries
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     outputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: rm(nlayer,nesp)   ! species volume mixing ratio 
-      
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer :: l, iq
-      
-c     tracer indexes in the chemistry:
-
-      integer,parameter :: i_co2  =  1
-      integer,parameter :: i_co   =  2
-      integer,parameter :: i_o    =  3
-      integer,parameter :: i_o1d  =  4
-      integer,parameter :: i_o2   =  5
-      integer,parameter :: i_o3   =  6
-      integer,parameter :: i_h    =  7
-      integer,parameter :: i_h2   =  8
-      integer,parameter :: i_oh   =  9
-      integer,parameter :: i_ho2  = 10
-      integer,parameter :: i_h2o2 = 11
-      integer,parameter :: i_ch4  = 12
-      integer,parameter :: i_h2o  = 13
-      integer,parameter :: i_n2   = 14
-      integer,parameter :: i_hox  = 15
-      integer,parameter :: i_ox   = 16
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialise chemical species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         rm(l,i_co2)  = max(zycol(l, igcm_co2),     1.e-30)
-         rm(l,i_co)   = max(zycol(l, igcm_co),      1.e-30)
-         rm(l,i_o)    = max(zycol(l, igcm_o),       1.e-30)
-         rm(l,i_o1d)  = max(zycol(l, igcm_o1d),     1.e-30)
-         rm(l,i_o2)   = max(zycol(l, igcm_o2),      1.e-30)
-         rm(l,i_o3)   = max(zycol(l, igcm_o3),      1.e-30)
-         rm(l,i_h)    = max(zycol(l, igcm_h),       1.e-30)
-         rm(l,i_h2)   = max(zycol(l, igcm_h2),      1.e-30)
-         rm(l,i_oh)   = max(zycol(l, igcm_oh),      1.e-30)
-         rm(l,i_ho2)  = max(zycol(l, igcm_ho2),     1.e-30)
-         rm(l,i_h2o2) = max(zycol(l, igcm_h2o2),    1.e-30)
-         rm(l,i_n2)   = max(zycol(l, igcm_n2),      1.e-30)
-         rm(l,i_h2o)  = max(zycol(l, igcm_h2o_vap), 1.e-30)
-      end do 
-
-      if (igcm_ch4 .eq. 0) then
-         do l = 1,lswitch-1
-            rm(l,i_ch4) = 0.
-         end do
-      else
-         do l = 1,lswitch-1
-            rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
-         end do
-      end if
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialise chemical families                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         rm(l,i_hox) = rm(l,i_h) 
-     $               + rm(l,i_oh)
-     $               + rm(l,i_ho2)
-         rm(l,i_ox)  = rm(l,i_o) 
-     $               + rm(l,i_o3)
-      end do
-
-      end subroutine gcmtochim
-
-c*****************************************************************
-c
-      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
-c
-c*****************************************************************
-c
-      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
-     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
-     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
-     &                      igcm_ch4
-
-      implicit none
-
-#include "callkeys.h"
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer, intent(in) :: nq     ! number of tracers
-      integer :: nesp               ! number of species in the chemistry
-      integer :: lswitch            ! interface level between chemistries
-      
-      real :: rm(nlayer,nesp)       ! species volume mixing ratio 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: zycol(nlayer,nq)  ! species volume mixing ratio in the gcm
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer l, iq
-      
-c     tracer indexes in the chemistry:
-
-      integer,parameter :: i_co2  =  1
-      integer,parameter :: i_co   =  2
-      integer,parameter :: i_o    =  3
-      integer,parameter :: i_o1d  =  4
-      integer,parameter :: i_o2   =  5
-      integer,parameter :: i_o3   =  6
-      integer,parameter :: i_h    =  7
-      integer,parameter :: i_h2   =  8
-      integer,parameter :: i_oh   =  9
-      integer,parameter :: i_ho2  = 10
-      integer,parameter :: i_h2o2 = 11
-      integer,parameter :: i_ch4  = 12
-      integer,parameter :: i_h2o  = 13
-      integer,parameter :: i_n2   = 14
-      integer,parameter :: i_hox  = 15
-      integer,parameter :: i_ox   = 16
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     save mixing ratios for the gcm
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         zycol(l, igcm_co2)     = rm(l,i_co2) 
-         zycol(l, igcm_co)      = rm(l,i_co) 
-         zycol(l, igcm_o)       = rm(l,i_o) 
-         zycol(l, igcm_o1d)     = rm(l,i_o1d)
-         zycol(l, igcm_o2)      = rm(l,i_o2) 
-         zycol(l, igcm_o3)      = rm(l,i_o3) 
-         zycol(l, igcm_h)       = rm(l,i_h)  
-         zycol(l, igcm_h2)      = rm(l,i_h2) 
-         zycol(l, igcm_oh)      = rm(l,i_oh) 
-         zycol(l, igcm_ho2)     = rm(l,i_ho2) 
-         zycol(l, igcm_h2o2)    = rm(l,i_h2o2)
-         zycol(l, igcm_n2)      = rm(l,i_n2)
-         zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
-      end do 
-
-      if (igcm_ch4 .ne. 0) then
-         do l = 1,lswitch-1
-            zycol(l,igcm_ch4) = rm(l,i_ch4)
-         end do
-      end if
-
-      end subroutine chimtogcm
-
-c*****************************************************************
-
-      subroutine chemrates(nlayer,
-     $                     lswitch, dens, press, t,
-     $                     surfdust1d, surfice1d,
-     $                     a001, a002, a003,
-     $                     b001, b002, b003, b004, b005, b006,
-     $                     b007, b008, b009,
-     $                     c001, c002, c003, c004, c005, c006,
-     $                     c007, c008, c009, c010, c011, c012,
-     $                     c013, c014, c015, c016, c017, c018,
-     $                     d001, d002, d003,
-     $                     e001, e002, e003,
-     $                     h001, h002, h003, h004, h005,
-     $                     t001, t002, t003, tau)
-
-c*****************************************************************
-
-      USE comcstfi_h
-      implicit none
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:                                                    c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer, intent(in) :: nlayer   ! number of atmospheric layers
-      integer :: lswitch              ! interface level between chemistries
-
-      real :: dens(nlayer)            ! density (cm-3)
-      real :: press(nlayer)           ! pressure (hpa)
-      real :: t(nlayer)               ! temperature (k)
-      real :: surfdust1d(nlayer)      ! dust surface area (cm^2/cm^3)
-      real :: surfice1d(nlayer)       ! ice surface area (cm^2/cm^3)
-      real :: tau                     ! dust opacity at 7 hpa
-
-      real, parameter :: tribo   = 0. ! switch for triboelectricity
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     outputs:                                                   c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: a001(nlayer), a002(nlayer), a003(nlayer)
-      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
-     $        b004(nlayer), b005(nlayer), b006(nlayer),
-     $        b007(nlayer), b008(nlayer), b009(nlayer)
-      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
-     $        c004(nlayer), c005(nlayer), c006(nlayer), 
-     $        c007(nlayer), c008(nlayer), c009(nlayer), 
-     $        c010(nlayer), c011(nlayer), c012(nlayer),
-     $        c013(nlayer), c014(nlayer), c015(nlayer), 
-     $        c016(nlayer), c017(nlayer), c018(nlayer)
-      real :: d001(nlayer), d002(nlayer), d003(nlayer)
-      real :: e001(nlayer), e002(nlayer), e003(nlayer)
-      real :: h001(nlayer), h002(nlayer), h003(nlayer),
-     $        h004(nlayer), h005(nlayer)
-      real :: t001(nlayer), t002(nlayer), t003(nlayer)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:                                                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
-      real :: ef, efmax, lossh2o, lossch4, lossco2
-
-      integer :: l
-      real :: k1a, k1b, k1a0, k1b0, k1ainf
-      real :: x, y, fc, fx
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     compute reaction rates
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oxygen compounds
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      a001: o + o2 + co2 -> o3 + co2
-c
-c        jpl 2003
-c
-c        co2 efficiency as a third body (2.075)
-c        from sehested et al., j. geophys. res., 100, 1995.
-c
-         a001(l) = 2.075
-     $             *6.0e-34*(t(l)/300.)**(-2.4)*dens(l)
-c
-c        mulcahy and williams, 1968
-c
-c        a001(l) = 2.68e-33*(t(l)/298.)**(-2.4)*dens(l)
-c
-c        nair et al., 1994
-c
-c        a001(l) = 1.3e-34*exp(724./t(l))*dens(l)
-c
-ccc      a002: o + o + co2 -> o2 + co2
-c
-c        Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
-c
-c        a002(l) = 2.5*5.2e-35*exp(900./t(l))*dens(l)
-c
-c        Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
-c
-c        a002(l) = 1.2e-32*(300./t(l))**(2.0)*dens(l)  ! yung expression
-c
-         a002(l) = 2.5*9.46e-34*exp(485./t(l))*dens(l) ! nist expression
-c
-c        baulch et al., 1976 confirmed by smith and robertson, 2008
-c
-c        a002(l) = 2.5*2.76e-34*exp(720./t(l))*dens(l)
-c
-ccc      a003: o + o3 -> o2 + o2
-c
-c        jpl 2003
-c
-         a003(l) = 8.0e-12*exp(-2060./t(l))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        reactions with o(1d)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      b001: o(1d) + co2  -> o + co2
-c
-c        jpl 2003
-c
-c        b001(l) = 7.4e-11*exp(120./t(l))
-c
-c        jpl 2006
-c
-         b001(l) = 7.5e-11*exp(115./t(l))
-c
-ccc      b002: o(1d) + h2o  -> oh + oh
-c
-c        jpl 2003
-c
-c        b002(l) = 2.2e-10
-c
-c        jpl 2006
-c 
-         b002(l) = 1.63e-10*exp(60./t(l))
-c    
-ccc      b003: o(1d) + h2  -> oh + h
-c
-c        jpl 2011
-c
-         b003(l) = 1.2e-10
-c    
-ccc      b004: o(1d) + o2  -> o + o2
-c
-c        jpl 2003
-c
-c        b004(l) = 3.2e-11*exp(70./t(l))
-c
-c        jpl 2006
-c
-         b004(l) = 3.3e-11*exp(55./t(l))
-c    
-ccc      b005: o(1d) + o3  -> o2 + o2
-c
-c        jpl 2003
-c
-         b005(l) = 1.2e-10
-c    
-ccc      b006: o(1d) + o3  -> o2 + o + o
-c
-c        jpl 2003
-c
-         b006(l) = 1.2e-10
-c    
-ccc      b007: o(1d) + ch4 -> ch3 + oh
-c
-c        jpl 2003
-c
-         b007(l) = 1.5e-10*0.75
-c
-ccc      b008: o(1d) + ch4 -> ch3o + h
-c
-c        jpl 2003
-c
-         b008(l) = 1.5e-10*0.20
-c
-ccc      b009: o(1d) + ch4 -> ch2o + h2
-c
-c        jpl 2003
-c
-         b009(l) = 1.5e-10*0.05
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        hydrogen compounds    
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      c001: o + ho2 -> oh + o2
-c
-c        jpl 2003
-c 
-         c001(l) = 3.0e-11*exp(200./t(l))
-c
-ccc      c002: o + oh -> o2 + h
-c
-c        jpl 2011
-c 
-         c002(l) = 1.8e-11*exp(180./t(l))
-c
-c        robertson and smith, j. chem. phys. a 110, 6673, 2006
-c
-c        c002(l) = 11.2e-11*t(l)**(-0.32)*exp(177./t(l))
-c
-ccc      c003: h + o3 -> oh + o2
-c
-c        jpl 2003
-c 
-         c003(l) = 1.4e-10*exp(-470./t(l))
-c
-ccc      c004: h + ho2 -> oh + oh
-c
-c        jpl 2003
-c 
-c        c004(l) = 8.1e-11*0.90
-c
-c        jpl 2006
-c
-         c004(l) = 7.2e-11
-c
-ccc      c005: h + ho2 -> h2 + o2
-c
-c        jpl 2003
-c 
-c        c005(l) = 8.1e-11*0.08
-c
-c        jpl 2006
-c
-         c005(l) = 6.9e-12
-c
-ccc      c006: h + ho2 -> h2o + o
-c
-c        jpl 2003
-c 
-c        c006(l) = 8.1e-11*0.02
-c
-c        jpl 2006
-c
-         c006(l) = 1.6e-12
-c
-ccc      c007: oh + ho2 -> h2o + o2
-c
-c        jpl 2003
-c 
-         c007(l) = 4.8e-11*exp(250./t(l))
-c
-c        jpl 2003 +20% d'apres canty et al., grl, 2006
-c 
-c        c007(l) = 4.8e-11*exp(250./t(l))*1.2
-c
-ccc      c008: ho2 + ho2 -> h2o2 + o2
-c
-c        jpl 2003
-c 
-c        c008(l) = 2.3e-13*exp(600./t(l))
-c
-c        christensen et al., grl, 13, 2002
-c
-         c008(l) = 1.5e-12*exp(19./t(l))
-c
-ccc      c009: oh + h2o2 -> h2o + ho2
-c
-c        jpl 2003
-c 
-c        c009(l) = 2.9e-12*exp(-160./t(l))
-c
-c        jpl 2006
-c
-         c009(l) = 1.8e-12
-c
-ccc      c010: oh + h2 -> h2o + h
-c
-c        jpl 2003
-c 
-c        c010(l) = 5.5e-12*exp(-2000./t(l))
-c
-c        jpl 2006
-c
-         c010(l) = 2.8e-12*exp(-1800./t(l))
-c
-ccc      c011: h + o2 + co2 -> ho2 + co2
-c
-c        jpl 2011
-c
-         ak0 = 2.5*4.4e-32*(t(l)/300.)**(-1.3)
-         ak1 = 7.5e-11*(t(l)/300.)**(0.2)
-c
-         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
-         c011(l) = rate*0.6**xpo
-c
-ccc      c012: o + h2o2 -> oh + ho2
-c
-c        jpl 2003
-c 
-         c012(l) = 1.4e-12*exp(-2000./t(l))
-c
-ccc      c013: oh + oh -> h2o + o
-c
-c        jpl 2003
-c 
-c        c013(l) = 4.2e-12*exp(-240./t(l))
-c
-c        jpl 2006
-c
-         c013(l) = 1.8e-12
-c
-ccc      c014: oh + o3 -> ho2 + o2
-c
-c        jpl 2003
-c
-         c014(l) = 1.7e-12*exp(-940./t(l))
-c
-c        jpl 2000
-c
-c        c014(l) = 1.5e-12*exp(-880./t(l))
-c
-c        nair et al., 1994 (jpl 1997)
-c
-c        c014(l) = 1.6e-12*exp(-940./t(l))
-c
-ccc      c015: ho2 + o3 -> oh + o2 + o2
-c
-c        jpl 2003
-c
-         c015(l) = 1.0e-14*exp(-490./t(l))
-c
-c        jpl 2000
-c
-c        c015(l) = 2.0e-14*exp(-680./t(l))
-c
-c        nair et al., 1994 (jpl 1997)
-c
-c        c015(l) = 1.1e-14*exp(-500./t(l))
-c
-ccc      c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
-c
-c        jpl 2011
-c
-         c016(l) = 2.5*2.1e-33
-     $              *exp(920./t(l))*dens(l)
-c
-ccc      c017: oh + oh + co2 -> h2o2 + co2
-c
-c        jpl 2003
-c
-         ak0 = 2.5*6.9e-31*(t(l)/300.)**(-1.0)
-         ak1 = 2.6e-11*(t(l)/300.)**(0.0)
-c
-c        jpl 1997
-c
-c        ak0 = 2.5*6.2e-31*(t(l)/300.)**(-1.0)
-c        ak1 = 2.6e-11*(t(l)/300.)**(0.0)
-c
-c        nair et al., 1994
-c
-c        ak0 = 2.5*7.1e-31*(t(l)/300.)**(-0.8)
-c        ak1 = 1.5e-11*(t(l)/300.)**(0.0)
-c
-         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
-         c017(l) = rate*0.6**xpo
-c
-ccc      c018: h + h + co2 -> h2 + co2
-c
-c        baulch et al., 1992
-c
-c        c018(l) = 2.5*8.85e-33*(t(l)/298.)**(-0.6)*dens(l)
-c
-c        baulch et al., 2005
-c
-         c018(l) = 2.5*1.8e-30*(t(l)**(-1.0))*dens(l)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        nitrogen compounds
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      d001: no2 + o -> no + o2
-c
-c        jpl 2003
-c
-c        d001(l) = 5.6e-12*exp(180./t(l))
-c
-ccc      jpl 2006
-c
-         d001(l) = 5.1e-12*exp(210./t(l))
-c
-ccc      d002: no + o3 -> no2 + o2
-c
-c        jpl 2003
-c
-         d002(l) = 3.0e-12*exp(-1500./t(l))
-c
-ccc      d003: no + ho2 -> no2 + oh
-c
-c        jpl 2011
-c
-         d003(l) = 3.3e-12*exp(270./t(l))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        carbon compounds
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      e001: oh + co -> co2 + h
-c
-c        jpl 2003
-c
-c        e001(l) = 1.5e-13*(1 + 0.6*press(l)/1013.)
-c
-c        mccabe et al., grl, 28, 3135, 2001
-c
-c        e001(l) = 1.57e-13 + 3.54e-33*dens(l)
-c
-c        atkinson et al. 2006
-c
-c        e001(l) = 1.44e-13 + 3.43e-33*dens(l)
-c
-c        joshi et al., 2006
-c
-         k1a0 = 1.34*2.5*dens(l)
-     $         *1/(1/(3.62e-26*t(l)**(-2.739)*exp(-20./t(l)))
-     $         + 1/(6.48e-33*t(l)**(0.14)*exp(-57./t(l))))   ! corrige de l'erreur publi
-         k1b0 = 1.17e-19*t(l)**(2.053)*exp(139./t(l))
-     $        + 9.56e-12*t(l)**(-0.664)*exp(-167./t(l))
-         k1ainf = 1.52e-17*t(l)**(1.858)*exp(28.8/t(l))
-     $          + 4.78e-8*t(l)**(-1.851)*exp(-318./t(l))
-         x = k1a0/(k1ainf - k1b0)
-         y = k1b0/(k1ainf - k1b0)
-         fc = 0.628*exp(-1223./t(l)) + (1. - 0.628)*exp(-39./t(l))
-     $      + exp(-t(l)/255.)
-         fx = fc**(1./(1. + (alog(x))**2))                   ! corrige de l'erreur publi
-         k1a = k1a0*((1. + y)/(1. + x))*fx
-         k1b = k1b0*(1./(1.+x))*fx
-c
-         e001(l) = k1a + k1b
-c
-ccc      e002: o + co + m -> co2 + m
-c
-c        tsang and hampson, 1986.
-c
-         e002(l) = 2.5*6.5e-33*exp(-2184./t(l))*dens(l)
-c
-c        baulch et al., butterworths, 1976.
-c
-c        e002(l) = 1.6e-32*exp(-2184./t(l))*dens(l)
-c
-ccc      e003: ch4 + oh -> ch3 + h2o
-c
-c        jpl 2003
-c
-c        e003(l) = 2.45e-12*exp(-1775./t(l))
-c
-c        jpl 2003, three-parameter expression
-c
-         e003(l) = 2.80e-14*(t(l)**0.667)*exp(-1575./t(l))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        heterogenous chemistry
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c        k = (surface*v*gamma)/4 (s-1)
-c        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
-c
-ccc      h001: ho2 + ice -> products
-c
-c        cooper and abbatt, 1996: gamma = 0.025
-c
-         h001(l) = surfice1d(l)
-     $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.025/4.
-c
-c        h002: oh + ice -> products
-c
-c        cooper and abbatt, 1996: gamma = 0.03
-c
-         h002(l) = surfice1d(l)
-     $            *100.*sqrt(8.*8.31*t(l)/(17.e-3*pi))*0.03/4.
-c
-c        h003: ho2 + dust -> products
-c
-c        jacob, 2000: gamma = 0.2
-c        see dereus et al., atm. chem. phys., 2005
-c
-c        h003(l) = surfdust1d(l)
-c    $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.2/4.
-         h003(l) = 0.     ! advised
-c
-ccc      h004: h2o2 + ice -> products
-c
-c        gamma = 1.e-3    test value
-c
-c        h004(l) = surfice1d(l)
-c    $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*0.001/4.
-         h004(l) = 0.     ! advised
-c
-c        h005: h2o2 + dust -> products
-c
-c        gamma = 5.e-4
-c        see dereus et al., atm. chem. phys., 2005
-c
-         h005(l) = surfdust1d(l)
-     $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*5.e-4/4.
-         h005(l) = 0.     ! advised
-c
-      end do
-c
-      if (tribo .eq. 1.) then
-c
-c        electrochemical reactions
-c
-c        efmax: maximum electric field (kv.m-1)
-c
-         efmax = 23.3
-c
-c        ef: actual electric field, scaled by tau.
-c
-c        if (tau .ge. 1.) then
-c           ef = efmax
-c        else
-c           ef = 0.
-c        end if
-c        ef = min(efmax,efmax*tau/1.0)
-c
-         ef = (efmax/0.5)*tau - (efmax/0.5)*0.5
-c
-         ef = max(ef, 0.)
-         ef = min(ef, efmax)
-c
-ccc      t001: h2o + e -> oh + h-
-c
-c        lossh2o: fit of oh/h- production rates
-c        given by delory et al., astrobiology, 6, 451, 2006
-c
-         if (ef .eq. 0.) then
-            lossh2o = 0.
-         else if (ef .lt. 10.) then
-            lossh2o = 0.054136*exp(1.0978*ef)
-         else if (ef .lt. 16.) then
-            lossh2o = 64.85*exp(0.38894*ef)
-         else if (ef .le. 20.) then
-            lossh2o = 0.2466*exp(0.73719*ef)
-         else
-            lossh2o = 2.3269e-8*exp(1.546*ef)
-         end if
-c
-c        production rates are given for h2o = 20 prec. microns.
-c        t001 is converted to first-order reaction rate
-c        assuming h2o number density at the surface =  5e13 mol cm-3
-c
-         do l = 1,21                     ! 70 km
-            t001(l) = lossh2o/5.e13      ! s-1
-         end do
-         do l = 22,lswitch-1
-            t001(l) = 0.
-         end do
-c
-ccc      t002: ch4 + e -> products
-c
-c        lossch4: fit of ch4 loss rates
-c        given by farrell et al., grl, 33, 2006
-c
-         if (ef .eq. 0.) then
-            lossch4 = 0.
-         else if (ef .gt. 20.) then
-            lossch4 = 1.113e-21*exp(1.6065*ef)
-         else if (ef .gt. 17.5) then
-            lossch4 = 1.e-15*exp(0.92103*ef)
-         else if (ef .gt. 14.) then
-            lossch4 = 1.e-13*exp(0.65788*ef)
-         else
-            lossch4 = 8.9238e-15*exp(0.835*ef)
-         end if
-c
-         do l = 1,21               ! 70 km
-            t002(l) = lossch4      ! s-1
-         end do
-         do l = 22,lswitch-1
-            t002(l) = 0.
-         end do
-c
-ccc      t003: co2 + e -> co + o-
-c
-c        lossco2: fit of co/o- production rates
-c        given by delory et al., astrobiology, 6, 451, 2006
-c
-         if (ef .eq. 0.) then
-            lossco2 = 0.
-         else if (ef .lt. 10.) then
-            lossco2 = 22.437*exp(1.045*ef)
-         else if (ef .lt. 16.) then
-            lossco2 = 17518.*exp(0.37896*ef)
-         else if (ef .lt. 20.) then
-            lossco2 = 54.765*exp(0.73946*ef)
-         else
-            lossco2 = 4.911e-6*exp(1.5508*ef)
-         end if
-c
-c        production rates are assumed to be given for p = 6 hPa
-c        lossco2 is converted to first-order reaction rate
-c        assuming co2 number density at the surface =  2e17 mol cm-3
-c
-         do l = 1,21                     ! 70 km
-            t003(l) = lossco2/2.e17      ! s-1
-         end do
-         do l = 22,lswitch-1
-            t003(l) = 0.
-         end do
-      else
-         do l = 1,lswitch-1
-            t001(l) = 0.
-            t002(l) = 0.
-            t003(l) = 0.
-         end do
-      end if
-c
-      end subroutine chemrates
-
-      end subroutine photochemistry
Index: /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2007)
+++ /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2007)
@@ -0,0 +1,1853 @@
+!*****************************************************************
+!
+!     Photochemical routine 
+!
+!     Author: Franck Lefevre
+!     ------
+!
+!     Version: 27/04/2017
+!
+!     ASIS scheme : for details on the method see
+!     Cariolle et al., Geosci. Model Dev., 10, 1467-1485, 2017.
+!
+!*****************************************************************
+
+subroutine photochemistry(nlayer, nq,                                 &
+                          ig, lswitch, zycol, sza, ptimestep, press,  &
+                          temp, dens, zmmean, dist_sol, surfdust1d,   &
+                          surfice1d, jo3, tau, iter)
+
+implicit none
+
+#include "chimiedata.h" 
+#include "callkeys.h"
+
+!===================================================================
+!     inputs:
+!===================================================================
+
+integer, intent(in) :: nlayer ! number of atmospheric layers
+integer, intent(in) :: nq     ! number of tracers
+integer :: ig                 ! grid point index
+      
+real :: sza                   ! solar zenith angle (deg)
+real :: ptimestep             ! physics timestep (s)
+real :: press(nlayer)         ! pressure (hpa)
+real :: temp(nlayer)          ! temperature (k)
+real :: dens(nlayer)          ! density (cm-3)
+real :: zmmean(nlayer)        ! mean molar mass (g/mole)
+real :: dist_sol              ! sun distance (au) 
+real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
+real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
+real :: tau                   ! optical depth at 7 hpa
+
+!===================================================================
+!     input/output:
+!===================================================================
+      
+real :: zycol(nlayer,nq)       ! chemical species volume mixing ratio
+
+!===================================================================
+!     output:
+!===================================================================
+      
+integer :: iter(nlayer)        ! iteration counter
+real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
+
+!===================================================================
+!     local:
+!===================================================================
+
+integer :: phychemrat         ! (physical timestep)/(nominal chemical timestep)
+integer :: j_o3_o1d, ilev
+integer :: iesp, nesp
+integer :: lswitch
+
+logical, save :: firstcall = .true.
+
+parameter (nesp = 17)         ! number of species in the chemistry
+
+! tracer indexes in the chemistry:
+
+integer,parameter :: i_co2  =  1
+integer,parameter :: i_co   =  2
+integer,parameter :: i_o    =  3
+integer,parameter :: i_o1d  =  4
+integer,parameter :: i_o2   =  5
+integer,parameter :: i_o3   =  6
+integer,parameter :: i_h    =  7
+integer,parameter :: i_h2   =  8
+integer,parameter :: i_oh   =  9
+integer,parameter :: i_ho2  = 10
+integer,parameter :: i_h2o2 = 11
+integer,parameter :: i_h2o  = 12
+integer,parameter :: i_n    = 13
+integer,parameter :: i_n2d  = 14
+integer,parameter :: i_no   = 15
+integer,parameter :: i_no2  = 16
+integer,parameter :: i_n2   = 17
+
+real :: ctimestep           ! standard timestep for the chemistry (s) 
+real :: dt_guess            ! first-guess timestep (s) 
+real :: dt_corrected        ! corrected timestep (s) 
+real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
+real :: time                ! internal time (between 0 and ptimestep, in s)
+
+real, dimension(nlayer,nesp)            :: rm   ! mixing ratios 
+real (kind = 8), dimension(nesp)        :: cold ! number densities at previous timestep (molecule.cm-3) 
+real (kind = 8), dimension(nlayer,nesp) :: c    ! number densities at current timestep (molecule.cm-3) 
+real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3) 
+ 
+! reaction rates
+ 
+real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
+real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
+real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
+logical :: hetero_dust, hetero_ice
+
+! matrix
+
+real (kind = 8), dimension(nesp,nesp) :: mat, mat1
+integer, dimension(nesp)              :: indx
+integer                               :: code
+
+! production and loss terms (for first-guess solution only)
+
+real (kind = 8), dimension(nesp) :: prod, loss
+
+!===================================================================
+!     initialisation of the reaction indexes
+!===================================================================
+
+if (firstcall) then
+   print*,'photochemistry: initialize indexes'
+   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
+               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
+               i_n, i_n2d, i_no, i_no2, i_n2)
+   firstcall = .false.
+end if
+
+!===================================================================
+!     initialisation of mixing ratios and densities       
+!===================================================================
+
+call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
+               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
+               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
+               i_n, i_n2d, i_no, i_no2, i_n2, dens, rm, c) 
+
+!===================================================================
+!     interpolation of photolysis rates in the lookup table      
+!===================================================================
+
+call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & 
+                rm(:,i_co2), rm(:,i_o3), v_phot)
+
+! save o3 photolysis for output
+
+j_o3_o1d = 5
+jo3(:) = v_phot(:,j_o3_o1d)
+
+!===================================================================
+!     reaction rates                                     
+!===================================================================
+!     switches for heterogeneous chemistry
+!     hetero_ice  : reactions on ice clouds
+!     hetero_dust : reactions on dust    
+!===================================================================
+
+hetero_dust = .false.
+hetero_ice  = .true.
+
+call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2), &
+                   press, temp, hetero_dust, hetero_ice,         &
+                   surfdust1d, surfice1d, v_phot, v_3, v_4)
+
+!===================================================================
+!     ctimestep : standard chemical timestep (s), defined as 
+!                 the fraction phychemrat of the physical timestep                           
+!===================================================================
+
+phychemrat = 1
+
+ctimestep = ptimestep/real(phychemrat)
+
+!print*, "ptimestep  = ", ptimestep
+!print*, "phychemrat = ", phychemrat
+!print*, "ctimestep  = ", ctimestep
+!stop
+
+!===================================================================
+!     loop over levels         
+!===================================================================
+
+do ilev = 1,lswitch - 1
+
+!  initializations
+
+   time = 0.
+   iter(ilev) = 0
+   dt_guess = ctimestep
+   cold(:) = c(ilev,:)
+
+!  internal loop for the chemistry
+
+   do while (time < ptimestep)
+
+   iter(ilev) = iter(ilev) + 1    ! iteration counter
+  
+!  first-guess: fill matrix
+
+   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
+
+!  adaptative evaluation of the sub time step
+
+   call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:),  &
+                  mat1, prod, loss, dens(ilev))
+
+   if (time + dt_corrected > ptimestep) then
+      dt_corrected = ptimestep - time
+   end if
+
+!  if (dt_corrected /= dt_guess) then  ! the timestep has been modified
+
+!  form the matrix identity + mat*dt_corrected
+
+   mat(:,:) = mat1(:,:)*dt_corrected
+   do iesp = 1,nesp
+      mat(iesp,iesp) = 1. + mat(iesp,iesp)
+   end do
+
+!  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
+     
+   cnew(:) = c(ilev,:)
+
+#ifdef LAPACK
+   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
+#else
+   write(*,*) "photochemistry error, missing LAPACK routine dgesv"
+   stop
+#endif
+
+!  end if
+
+!  eliminate small values
+
+   where (cnew(:)/dens(ilev) < 1.e-30)
+      cnew(:) = 0.
+   end where
+
+!  update concentrations
+
+   cold(:)   = c(ilev,:)
+   c(ilev,:) = cnew(:)
+   cnew(:)   = 0.
+
+!  increment internal time
+
+   time = time + dt_corrected
+   dt_guess = dt_corrected     ! first-guess timestep for next iteration
+
+   end do ! while (time < ptimestep)
+
+end do ! ilev
+
+!===================================================================
+!     save chemical species for the gcm       
+!===================================================================
+
+call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
+               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
+               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
+               i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
+
+contains
+
+!================================================================
+
+ subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
+                      prod, loss, dens)
+
+!================================================================
+! iterative evaluation of the appropriate time step dtnew
+! according to curvature criterion based on
+! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
+! with r = (tn - tn-1)/(tn+1 - tn)
+!================================================================
+
+implicit none
+
+! input
+
+integer :: nesp  ! number of species in the chemistry
+
+real :: dtold, ctimestep
+real (kind = 8), dimension(nesp)      :: cold, ccur
+real (kind = 8), dimension(nesp,nesp) :: mat1
+real (kind = 8), dimension(nesp)      :: prod, loss
+real                                  :: dens
+
+! output
+
+real :: dtnew
+
+! local
+
+real (kind = 8), dimension(nesp)      :: cnew
+real (kind = 8), dimension(nesp,nesp) :: mat
+real (kind = 8) :: atol, ratio, e, es, coef
+
+integer                  :: code, iesp, iter
+integer, dimension(nesp) :: indx
+
+real :: dttest
+
+! parameters
+
+real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
+real (kind = 8), parameter :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
+real (kind = 8), parameter :: rtol    = 1./0.05   ! 1/rtol recommended value : 0.1-0.02
+integer,         parameter :: niter   = 3        ! number of iterations
+real (kind = 8), parameter :: coefmax = 2.
+real (kind = 8), parameter :: coefmin = 0.1 
+logical                    :: fast_guess = .true.
+
+dttest = dtold   ! dttest = dtold = dt_guess
+
+atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
+
+do iter = 1,niter
+
+if (fast_guess) then
+
+! first guess : fast semi-implicit method
+
+   do iesp = 1, nesp
+      cnew(iesp) = (ccur(iesp) + prod(iesp)*dttest)/(1. + loss(iesp)*dttest)
+   end do
+
+else
+
+! first guess : form the matrix identity + mat*dt_guess
+
+   mat(:,:) = mat1(:,:)*dttest
+   do iesp = 1,nesp
+      mat(iesp,iesp) = 1. + mat(iesp,iesp)
+   end do
+
+! form right-hand side (RHS) of the system
+
+   cnew(:) = ccur(:)
+
+! solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
+
+#ifdef LAPACK
+      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
+#else
+   write(*,*) "photochemistry error, missing LAPACK routine dgesv"
+   stop
+#endif
+
+end if
+
+! ratio old/new subtimestep
+
+ratio = dtold/dttest
+
+! e : local error indicatocitr
+
+e = 0.
+
+do iesp = 1,nesp
+   es = 2.*abs((ratio*cnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp))   &
+         /(1. + ratio)/max(ccur(iesp),atol))
+
+   if (es > e) then
+      e = es
+   end if
+end do
+e = rtol*e
+
+! timestep correction
+
+coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
+
+dttest = max(dtmin,dttest*coef)
+dttest = min(ctimestep,dttest)
+
+end do ! iter
+
+! new timestep
+
+dtnew = dttest
+
+end subroutine define_dt
+
+!======================================================================
+
+ subroutine reactionrates(nlayer,                               &
+                          lswitch, dens, co2, o2, press, t,     &
+                          hetero_dust, hetero_ice,              &
+                          surfdust1d, surfice1d,                &
+                          v_phot, v_3, v_4)
+ 
+!================================================================
+! compute reaction rates                                        !
+!----------------------------------------------------------------
+! reaction               type                array              !
+!----------------------------------------------------------------
+! A + B    --> C + D     bimolecular         v_4                !
+! A + A    --> B + C     quadratic           v_3                !
+! A + C    --> B + C     quenching           v_phot             !
+! A + ice  --> B + C     heterogeneous       v_phot             !
+!================================================================
+
+use comcstfi_h
+
+implicit none
+
+#include "chimiedata.h"
+
+!----------------------------------------------------------------------
+!     input
+!----------------------------------------------------------------------
+
+integer, intent(in)     :: nlayer            ! number of atmospheric layers
+integer                 :: lswitch           ! interface level between lower
+                                             ! atmosphere and thermosphere chemistries
+real, dimension(nlayer) :: dens              ! total number density (molecule.cm-3)
+real, dimension(nlayer) :: press             ! pressure (hPa)
+real, dimension(nlayer) :: t                 ! temperature (K)
+real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
+real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
+real (kind = 8), dimension(nlayer) :: co2    ! co2 number density (molecule.cm-3)
+real (kind = 8), dimension(nlayer) :: o2     ! o2 number density (molecule.cm-3)
+logical :: hetero_dust, hetero_ice           ! switches for heterogeneous chemistry
+
+!----------------------------------------------------------------------
+!     output
+!----------------------------------------------------------------------
+
+real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
+real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
+real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
+
+!----------------------------------------------------------------------
+!     local
+!----------------------------------------------------------------------
+
+integer          :: ilev
+integer          :: nb_phot, nb_reaction_3, nb_reaction_4
+real :: ak0, ak1, xpo, rate
+real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
+real, dimension(nlayer) :: deq
+real, dimension(nlayer) :: a001, a002, a003,                           &
+                             b001, b002, b003, b004, b005, b006, b007,   &
+                             b008, b009,                                 &
+                             c001, c002, c003, c004, c005, c006, c007,   &
+                             c008, c009, c010, c011, c012, c013, c014,   &
+                             c015, c016, c017, c018,                     &
+                             d001, d002, d003, d004, d005, d006, d007,   &
+                             d008, d009,                                 &
+                             e001, e002,                                 &
+                             h001, h002, h003, h004, h005
+
+!----------------------------------------------------------------------
+!     initialisation
+!----------------------------------------------------------------------
+
+      nb_phot       = 11       ! jmars.20140930 reduit de 13 a 11
+      nb_reaction_3 = 0
+      nb_reaction_4 = 0
+
+!----------------------------------------------------------------------
+!     oxygen reactions 
+!----------------------------------------------------------------------
+
+!---  a001: o + o2 + co2 -> o3 + co2
+
+!     jpl 2003
+!
+!     co2/n2 efficiency as a third body = 2.075
+!     from sehested et al., j. geophys. res., 100, 1995.
+
+      a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = a001(:)
+
+!---  a002: o + o + co2 -> o2 + co2
+
+!     Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
+
+!     a002(:) = 2.5*5.2e-35*exp(900./t(:))*dens(:)
+
+!     Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
+
+!     a002(:) = 1.2e-32*(300./t(:))**(2.0)*dens(:)  ! yung expression
+      a002(:) = 2.5*9.46e-34*exp(485./t(:))*dens(:) ! nist expression
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = a002(:)
+
+!---  a003: o + o3 -> o2 + o2
+
+!     jpl 2003
+
+      a003(:) = 8.0e-12*exp(-2060./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = a003(:)
+
+!----------------------------------------------------------------------
+!     o(1d) reactions
+!----------------------------------------------------------------------
+
+!---  b001: o(1d) + co2  -> o + co2
+
+!     jpl 2006
+
+      b001(:) = 7.5e-11*exp(115./t(:))
+   
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = b001(:)*co2(:)
+
+!---  b002: o(1d) + h2o  -> oh + oh
+
+!     jpl 2006
+ 
+      b002(:) = 1.63e-10*exp(60./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = b002(:)
+
+!---  b003: o(1d) + h2  -> oh + h
+
+!     jpl 2011
+
+      b003(:) = 1.2e-10
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = b003(:)
+
+!---  b004: o(1d) + o2  -> o + o2
+
+!     jpl 2006
+
+      b004(:) = 3.3e-11*exp(55./t(:))
+
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = b004(:)*o2(:)
+    
+!---  b005: o(1d) + o3  -> o2 + o2
+
+!     jpl 2003
+
+      b005(:) = 1.2e-10
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = b005(:)
+    
+!---  b006: o(1d) + o3  -> o2 + o + o
+
+!     jpl 2003
+
+      b006(:) = 1.2e-10
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = b006(:)
+    
+!---  b007: o(1d) + ch4 -> ch3 + oh
+
+!     jpl 2003
+
+      b007(:) = 1.5e-10*0.75
+
+!---  b008: o(1d) + ch4 -> ch3o + h
+
+!     jpl 2003
+
+      b008(:) = 1.5e-10*0.20
+!
+!---  b009: o(1d) + ch4 -> ch2o + h2
+
+!     jpl 2003
+
+      b009(:) = 1.5e-10*0.05
+
+!----------------------------------------------------------------------
+!     hydrogen reactions
+!----------------------------------------------------------------------
+
+!---  c001: o + ho2 -> oh + o2
+
+!     jpl 2003
+
+      c001(:) = 3.0e-11*exp(200./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c001(:)
+
+!---  c002: o + oh -> o2 + h
+
+!     jpl 2011
+
+      c002(:) = 1.8e-11*exp(180./t(:))
+
+!     robertson and smith, j. chem. phys. a 110, 6673, 2006
+
+!     c002(:) = 11.2e-11*t(:)**(-0.32)*exp(177./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c002(:)
+
+!---  c003: h + o3 -> oh + o2
+
+!     jpl 2003
+
+      c003(:) = 1.4e-10*exp(-470./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c003(:)
+
+!---  c004: h + ho2 -> oh + oh
+
+!     jpl 2006
+
+      c004(:) = 7.2e-11
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c004(:)
+
+!---  c005: h + ho2 -> h2 + o2
+
+!     jpl 2006
+
+      c005(:) = 6.9e-12
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c005(:)
+
+!---  c006: h + ho2 -> h2o + o
+
+!     jpl 2006
+
+      c006(:) = 1.6e-12
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c006(:)
+
+!---  c007: oh + ho2 -> h2o + o2
+
+!     jpl 2003
+
+!     canty et al., grl, 2006 suggest to increase this rate
+!     by 20%. not done here.
+
+      c007(:) = 4.8e-11*exp(250./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c007(:)
+
+!---  c008: ho2 + ho2 -> h2o2 + o2
+
+!     jpl 2006
+
+!     c008(:) = 3.5e-13*exp(430./t(:))
+
+!     christensen et al., grl, 13, 2002
+
+      c008(:) = 1.5e-12*exp(19./t(:))
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = c008(:)
+
+!---  c009: oh + h2o2 -> h2o + ho2
+
+!     jpl 2006
+
+      c009(:) = 1.8e-12
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c009(:)
+
+!---  c010: oh + h2 -> h2o + h
+
+!     jpl 2006
+
+      c010(:) = 2.8e-12*exp(-1800./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c010(:)
+
+!---  c011: h + o2 + co2 -> ho2 + co2
+
+!     jpl 2011
+!     co2/n2 efficiency as a third body = 2.4
+!     from ashman and haynes, 27th symposium on combustion, 1998.
+
+      do ilev = 1,lswitch-1
+         ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3)
+         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
+
+         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
+         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
+         c011(ilev) = rate*0.6**xpo
+      end do
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c011(:)
+
+!---  c012: o + h2o2 -> oh + ho2
+
+!     jpl 2003
+
+      c012(:) = 1.4e-12*exp(-2000./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c012(:)
+
+!---  c013: oh + oh -> h2o + o
+
+!     jpl 2006
+
+      c013(:) = 1.8e-12
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = c013(:)
+
+!---  c014: oh + o3 -> ho2 + o2
+
+!     jpl 2003
+
+      c014(:) = 1.7e-12*exp(-940./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c014(:)
+
+!---  c015: ho2 + o3 -> oh + o2 + o2
+
+!     jpl 2003
+
+      c015(:) = 1.0e-14*exp(-490./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = c015(:)
+
+!---  c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
+
+!     jpl 2011
+
+      c016(:) = 2.5*2.1e-33*exp(920./t(:))*dens(:)
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = c016(:)
+
+!---  c017: oh + oh + co2 -> h2o2 + co2
+
+!     jpl 2003
+
+      do ilev = 1,lswitch-1
+         ak0 = 2.5*6.9e-31*(t(ilev)/300.)**(-1.0)
+         ak1 = 2.6e-11*(t(ilev)/300.)**(0.0)
+
+         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
+         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
+         c017(ilev) = rate*0.6**xpo
+      end do
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = c017(:)
+
+!---  c018: h + h + co2 -> h2 + co2
+
+!     baulch et al., 2005
+
+      c018(:) = 2.5*1.8e-30*(t(:)**(-1.0))*dens(:)
+
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = c018(:)
+
+!----------------------------------------------------------------------
+!     nitrogen reactions
+!----------------------------------------------------------------------
+
+!---  d001: no2 + o -> no + o2
+
+!     jpl 2006
+
+      d001(:) = 5.1e-12*exp(210./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d001(:)
+
+!---  d002: no + o3 -> no2 + o2
+
+!     jpl 2006
+
+      d002(:) = 3.0e-12*exp(-1500./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d002(:)
+
+!---  d003: no + ho2 -> no2 + oh
+
+!     jpl 2011
+
+      d003(:) = 3.3e-12*exp(270./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d003(:)
+
+!---  d004: n + no -> n2 + o
+
+!     jpl 2011
+
+      d004(:) = 2.1e-11*exp(100./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d004(:)
+
+!---  d005: n + o2 -> no + o
+
+!     jpl 2011
+
+      d005(:) = 1.5e-11*exp(-3600./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d005(:)
+
+!---  d006: no2 + h -> no + oh
+
+!     jpl 2011
+
+      d006(:) = 4.0e-10*exp(-340./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d006(:)
+
+!---  d007: n + o -> no
+ 
+      d007(:) = 2.8e-17*(300./t(:))**0.5
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d007(:)
+
+!---  d008: n + ho2 -> no + oh
+
+!     brune et al., j. chem. phys., 87, 1983 
+
+      d008(:) = 2.19e-11
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d008(:)
+
+!---  d009: n + oh -> no + h
+
+!     atkinson et al., j. phys. chem. ref. data, 18, 881, 1989 
+
+      d009(:) = 3.8e-11*exp(85./t(:))
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = d009(:)
+
+!----------------------------------------------------------------------
+!     carbon reactions
+!----------------------------------------------------------------------
+
+!---  e001: oh + co -> co2 + h
+
+!     jpl 2003
+
+!     e001(:) = 1.5e-13*(1 + 0.6*press(:)/1013.)
+
+!     mccabe et al., grl, 28, 3135, 2001
+
+!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
+
+!     jpl 2006
+
+!     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
+!     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
+!     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
+!     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
+
+!     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
+!     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
+!     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
+!     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
+
+!     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
+
+!     joshi et al., 2006
+
+      do ilev = 1,lswitch-1
+         k1a0 = 1.34*2.5*dens(ilev)                                  &
+               *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
+               + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
+         k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
+              + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
+         k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
+                + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
+         x = k1a0/(k1ainf - k1b0)
+         y = k1b0/(k1ainf - k1b0)
+         fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
+            + exp(-t(ilev)/255.)
+         fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
+         k1a = k1a0*((1. + y)/(1. + x))*fx
+         k1b = k1b0*(1./(1.+x))*fx
+
+         e001(ilev) = k1a + k1b
+      end do
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = e001(:)
+
+!---  e002: o + co + m -> co2 + m
+
+!     tsang and hampson, 1986.
+
+      e002(:) = 2.5*6.5e-33*exp(-2184./t(:))*dens(:)
+
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = e002(:)
+
+!----------------------------------------------------------------------
+!     heterogeneous chemistry 
+!----------------------------------------------------------------------
+
+      if (hetero_ice) then
+
+!        k = (surface*v*gamma)/4 (s-1)
+!        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
+ 
+!---     h001: ho2 + ice -> products
+ 
+!        cooper and abbatt, 1996: gamma = 0.025
+      
+         gam = 0.025
+         h001(:) = surfice1d(:)       &
+                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
+ 
+!        h002: oh + ice -> products
+ 
+!        cooper and abbatt, 1996: gamma = 0.03
+ 
+         gam = 0.03
+         h002(:) = surfice1d(:)       &
+                   *100.*sqrt(8.*8.31*t(:)/(17.e-3*pi))*gam/4.
+
+!---     h003: h2o2 + ice -> products
+ 
+!        gamma = 0.    test value
+ 
+         gam = 0.
+         h003(:) = surfice1d(:)       &
+                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
+      else
+         h001(:) = 0.
+         h002(:) = 0.
+         h003(:) = 0.
+      end if
+
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = h001(:)
+
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = h002(:)
+
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = h003(:)
+
+      if (hetero_dust) then
+ 
+!---     h004: ho2 + dust -> products
+ 
+!        jacob, 2000: gamma = 0.2
+!        see dereus et al., atm. chem. phys., 2005
+ 
+         gam = 0.2
+         h004(:) = surfdust1d(:)  &
+                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
+ 
+!---     h005: h2o2 + dust -> products
+ 
+!        gamma = 5.e-4
+!        see dereus et al., atm. chem. phys., 2005
+ 
+         gam = 5.e-4
+         h005(:) = surfdust1d(:)  &
+                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
+      else
+         h004(:) = 0.
+         h005(:) = 0.
+      end if
+ 
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = h004(:)
+
+      nb_phot = nb_phot + 1
+      v_phot(:,nb_phot) = h005(:)
+
+end subroutine reactionrates
+
+!======================================================================
+
+ subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
+
+!======================================================================
+! filling of the jacobian matrix
+!======================================================================
+
+use types_asis
+
+implicit none
+
+#include "chimiedata.h"
+
+! input
+
+integer             :: ilev    ! level index
+integer             :: nesp    ! number of species in the chemistry
+integer, intent(in) :: nlayer  ! number of atmospheric layers
+
+real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
+real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
+real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
+real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
+
+! output
+
+real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
+real (kind = 8), dimension(nesp), intent(out)      :: prod, loss
+
+! local
+
+integer :: iesp
+integer :: ind_phot_2,ind_phot_4,ind_phot_6
+integer :: ind_3_2,ind_3_4,ind_3_6
+integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
+integer :: iphot,i3,i4
+
+real(kind = 8) :: eps, eps_4  ! implicit/explicit coefficient
+
+! initialisations 
+
+mat(:,:) = 0.
+prod(:)  = 0.
+loss(:)  = 0.
+
+! photodissociations
+! or reactions a + c -> b + c
+! or reactions a + ice -> b + c
+
+do iphot = 1,nb_phot_max
+
+  ind_phot_2 = indice_phot(iphot)%z2
+  ind_phot_4 = indice_phot(iphot)%z4
+  ind_phot_6 = indice_phot(iphot)%z6
+
+  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
+  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
+  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
+
+  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
+  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+
+end do
+
+! reactions a + a -> b + c 
+
+do i3 = 1,nb_reaction_3_max
+
+  ind_3_2 = indice_3(i3)%z2
+  ind_3_4 = indice_3(i3)%z4
+  ind_3_6 = indice_3(i3)%z6
+
+  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
+  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)
+  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)
+
+  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
+  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+
+end do
+
+! reactions a + b -> c + d 
+
+eps = 1.d-10
+
+do i4 = 1,nb_reaction_4_max
+
+  ind_4_2 = indice_4(i4)%z2
+  ind_4_4 = indice_4(i4)%z4
+  ind_4_6 = indice_4(i4)%z6
+  ind_4_8 = indice_4(i4)%z8
+
+  eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
+  eps_4 = min(eps_4,1.0)
+
+  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4) 
+  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
+  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)   
+  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
+  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
+  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)
+  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+
+end do
+
+end subroutine fill_matrix
+
+!================================================================
+
+ 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)
+
+!================================================================
+! set the "indice" arrays used to fill the jacobian matrix      !
+!----------------------------------------------------------------
+! reaction                                   type               !
+!----------------------------------------------------------------
+! A + hv   --> B + C     photolysis          indice_phot        ! 
+! A + B    --> C + D     bimolecular         indice_4           !
+! A + A    --> B + C     quadratic           indice_3           !
+! A + C    --> B + C     quenching           indice_phot        !
+! A + ice  --> B + C     heterogeneous       indice_phot        !
+!================================================================
+
+use types_asis
+
+implicit none
+
+#include "chimiedata.h"
+
+! input
+
+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
+
+! local
+
+integer :: nb_phot, nb_reaction_3, nb_reaction_4
+integer :: i_dummy
+
+allocate (indice_phot(nb_phot_max))
+allocate (indice_3(nb_reaction_3_max))
+allocate (indice_4(nb_reaction_4_max))
+
+i_dummy = 1
+
+nb_phot       = 0
+nb_reaction_3 = 0
+nb_reaction_4 = 0
+
+!===========================================================
+!      O2 + hv -> O + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy) 
+
+!===========================================================
+!      O2 + hv -> O + O(1D)
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d) 
+
+!===========================================================
+!      CO2 + hv -> CO + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o) 
+
+!===========================================================
+!      CO2 + hv -> CO + O(1D)
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d) 
+
+!===========================================================
+!      O3 + hv -> O2 + O(1D)
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d) 
+
+!===========================================================
+!      O3 + hv -> O2 + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o) 
+
+!===========================================================
+!      H2O + hv -> H + OH
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh) 
+
+!===========================================================
+!      H2O2 + hv -> OH + OH
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy) 
+
+!===========================================================
+!      HO2 + hv -> OH + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o) 
+
+!===========================================================
+!      NO + hv -> N + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o) 
+
+!===========================================================
+!      NO2 + hv -> NO + O
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o) 
+
+!===========================================================
+!      a001 : O + O2 + CO2 -> O3 + CO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy) 
+
+!===========================================================
+!      a002 : O + O + CO2 -> O2 + CO2
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      a003 : O + O3 -> O2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      b001 : O(1D) + CO2 -> O + CO2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
+
+!===========================================================
+!      b002 : O(1D) + H2O -> OH + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy) 
+
+!===========================================================
+!      b003 : O(1D) + H2 -> OH + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h) 
+
+!===========================================================
+!      b004 : O(1D) + O2 -> O + O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
+
+!===========================================================
+!      b005 : O(1D) + O3 -> O2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
+
+!===========================================================
+!      b006 : O(1D) + O3 -> O2 + O + O
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 
+
+!===========================================================
+!      c001 : O + HO2 -> OH + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2) 
+
+!===========================================================
+!      c002 : O + OH -> O2 + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h) 
+
+!===========================================================
+!      c003 : H + O3 -> OH + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2) 
+
+!===========================================================
+!      c004 : H + HO2 -> OH + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy) 
+
+!===========================================================
+!      c005 : H + HO2 -> H2 + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2) 
+
+!===========================================================
+!      c006 : H + HO2 -> H2O + O
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o) 
+
+!===========================================================
+!      c007 : OH + HO2 -> H2O + O2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2) 
+
+!===========================================================
+!      c008 : HO2 + HO2 -> H2O2 + O2
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
+
+!===========================================================
+!      c009 : OH + H2O2 -> H2O + HO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2) 
+
+!===========================================================
+!      c010 : OH + H2 -> H2O + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h) 
+
+!===========================================================
+!      c011 : H + O2 + CO2 -> HO2 + CO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy) 
+
+!===========================================================
+!      c012 : O + H2O2 -> OH + HO2
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2) 
+
+!===========================================================
+!      c013 : OH + OH -> H2O + O
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o) 
+
+!===========================================================
+!      c014 : OH + O3 -> HO2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2) 
+
+!===========================================================
+!      c015 : HO2 + O3 -> OH + O2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2) 
+
+!===========================================================
+!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
+
+!===========================================================
+!      c017 : OH + OH + CO2 -> H2O2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy) 
+
+!===========================================================
+!      c018 : H + H + CO2 -> H2 + CO2 
+!===========================================================
+
+nb_reaction_3 = nb_reaction_3 + 1
+
+indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 
+
+!===========================================================
+!      d001 : NO2 + O -> NO + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2) 
+
+!===========================================================
+!      d002 : NO + O3 -> NO2 + O2 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2) 
+
+!===========================================================
+!      d003 : NO + HO2 -> NO2 + OH 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh) 
+
+!===========================================================
+!      d004 : N + NO -> N2 + O 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o) 
+
+!===========================================================
+!      d005 : N + O2 -> NO + O 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o) 
+
+!===========================================================
+!      d006 : NO2 + H -> NO + OH 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh) 
+
+!===========================================================
+!      d007 : N + O -> NO
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy) 
+
+!===========================================================
+!      d008 : N + HO2 -> NO + OH
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh) 
+
+!===========================================================
+!      d009 : N + OH -> NO + H
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h) 
+
+!===========================================================
+!      e001 : CO + OH -> CO2 + H 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h) 
+
+!===========================================================
+!      e002 : CO + O + M -> CO2 + M 
+!===========================================================
+
+nb_reaction_4 = nb_reaction_4 + 1
+
+indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 
+
+!===========================================================
+!      h001: HO2 + ice -> products
+!            treated as
+!            HO2 -> 0.5 H2O + 0.75 O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
+
+!===========================================================
+!      h002: OH + ice -> products
+!            treated as
+!            OH -> 0.5 H2O + 0.25 O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2) 
+
+!===========================================================
+!      h003: H2O2 + ice -> products
+!            treated as
+!            H2O2 -> H2O + 0.5 O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
+
+!===========================================================
+!      h004: HO2 + dust -> products
+!            treated as
+!            HO2 -> 0.5 H2O + 0.75 O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
+
+!===========================================================
+!      h005: H2O2 + dust -> products
+!            treated as
+!            H2O2 -> H2O + 0.5 O2
+!===========================================================
+
+nb_phot = nb_phot + 1
+
+indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
+
+!===========================================================
+!  check dimensions 
+!===========================================================
+
+print*, 'nb_phot       = ', nb_phot
+print*, 'nb_reaction_4 = ', nb_reaction_4
+print*, 'nb_reaction_3 = ', nb_reaction_3
+
+if ((nb_phot /= nb_phot_max)             .or.  &
+    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
+    (nb_reaction_4 /= nb_reaction_4_max)) then
+   print*, 'wrong dimensions in indice' 
+   stop
+end if  
+
+end subroutine indice
+
+!*****************************************************************
+
+      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
+                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
+                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
+                           i_n, i_n2d, i_no, i_no2, i_n2,            &
+                           dens, rm, c) 
+
+!*****************************************************************
+
+      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
+     &                      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
+
+      implicit none
+
+#include "callkeys.h"
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     input:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer, intent(in) :: nq     ! number of tracers in the gcm
+      integer :: nesp               ! number of species in the chemistry
+      integer :: lswitch            ! interface level between chemistries
+
+      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
+
+      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
+      real :: dens(nlayer)          ! total number density (molecule.cm-3) 
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     output:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      real, dimension(nlayer,nesp)            :: rm ! volume mixing ratios 
+      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 
+      
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     local:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      integer      :: l, iesp
+      logical,save :: firstcall = .true.
+      
+      
+!     first call initializations
+
+      if (firstcall) then
+
+!       identify the indexes of the tracers we need
+
+         if (igcm_co2 == 0) then
+            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
+            stop
+         endif
+         if (igcm_co == 0) then
+            write(*,*) "gcmtochim: Error; no CO tracer !!!"
+            stop
+         end if
+         if (igcm_o == 0) then
+            write(*,*) "gcmtochim: Error; no O tracer !!!"
+            stop
+         end if
+         if (igcm_o1d == 0) then
+            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
+            stop
+         end if
+         if (igcm_o2 == 0) then
+            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
+            stop
+         end if
+         if (igcm_o3 == 0) then
+            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
+            stop
+         end if
+         if (igcm_h == 0) then
+            write(*,*) "gcmtochim: Error; no H tracer !!!"
+            stop
+         end if
+         if (igcm_h2 == 0) then
+            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
+            stop
+         end if
+         if (igcm_oh == 0) then
+            write(*,*) "gcmtochim: Error; no OH tracer !!!"
+            stop
+         end if
+         if (igcm_ho2 == 0) then
+            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
+            stop
+         end if
+         if (igcm_h2o2 == 0) then
+            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
+            stop
+         end if
+         if (igcm_n == 0) then
+            write(*,*) "gcmtochim: Error; no N tracer !!!"
+            stop
+         end if
+         if (igcm_n2d == 0) then
+            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
+            stop
+         end if
+         if (igcm_no == 0) then
+            write(*,*) "gcmtochim: Error; no NO tracer !!!"
+            stop
+         end if
+         if (igcm_no2 == 0) then
+            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
+            stop
+         end if
+         if (igcm_n2 == 0) then
+            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
+            stop
+         end if
+         if (igcm_h2o_vap == 0) then
+            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
+            stop
+         end if
+         firstcall = .false.
+      end if ! of if (firstcall)
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     initialise mixing ratios
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      do l = 1,lswitch-1
+         rm(l,i_co2)  = zycol(l, igcm_co2)
+         rm(l,i_co)   = zycol(l, igcm_co)
+         rm(l,i_o)    = zycol(l, igcm_o)
+         rm(l,i_o1d)  = zycol(l, igcm_o1d)
+         rm(l,i_o2)   = zycol(l, igcm_o2)
+         rm(l,i_o3)   = zycol(l, igcm_o3)
+         rm(l,i_h)    = zycol(l, igcm_h)
+         rm(l,i_h2)   = zycol(l, igcm_h2)
+         rm(l,i_oh)   = zycol(l, igcm_oh)
+         rm(l,i_ho2)  = zycol(l, igcm_ho2)
+         rm(l,i_h2o2) = zycol(l, igcm_h2o2)
+         rm(l,i_h2o)  = zycol(l, igcm_h2o_vap)
+         rm(l,i_n)    = zycol(l, igcm_n)
+         rm(l,i_n2d)  = zycol(l, igcm_n2d)
+         rm(l,i_no)   = zycol(l, igcm_no)
+         rm(l,i_no2)  = zycol(l, igcm_no2)
+         rm(l,i_n2)   = zycol(l, igcm_n2)
+      end do 
+
+      where (rm(:,:) < 1.e-30)
+         rm(:,:) = 0.
+      end where
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     initialise number densities
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      do iesp = 1,nesp
+         do l = 1,lswitch-1
+            c(l,iesp) = rm(l,iesp)*dens(l)
+         end do
+      end do
+
+      end subroutine gcmtochim
+
+!*****************************************************************
+ 
+      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,         &
+                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
+                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
+                           i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
+ 
+!*****************************************************************
+ 
+      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,          &
+                            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
+
+      implicit none
+
+#include "callkeys.h"
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     inputs:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 
+      integer, intent(in) :: nlayer  ! number of atmospheric layers
+      integer, intent(in) :: nq      ! number of tracers in the gcm
+      integer :: nesp                ! number of species in the chemistry
+      integer :: lswitch             ! interface level between chemistries
+      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
+
+      real :: dens(nlayer)     ! total number density (molecule.cm-3) 
+      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     output:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+       
+      real zycol(nlayer,nq)  ! volume mixing ratios in the gcm
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     local:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      integer l
+      
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     save mixing ratios for the gcm
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      do l = 1,lswitch-1
+         zycol(l, igcm_co2)     = c(l,i_co2)/dens(l) 
+         zycol(l, igcm_co)      = c(l,i_co)/dens(l) 
+         zycol(l, igcm_o)       = c(l,i_o)/dens(l) 
+         zycol(l, igcm_o1d)     = c(l,i_o1d)/dens(l)
+         zycol(l, igcm_o2)      = c(l,i_o2)/dens(l) 
+         zycol(l, igcm_o3)      = c(l,i_o3)/dens(l) 
+         zycol(l, igcm_h)       = c(l,i_h)/dens(l)  
+         zycol(l, igcm_h2)      = c(l,i_h2)/dens(l) 
+         zycol(l, igcm_oh)      = c(l,i_oh)/dens(l) 
+         zycol(l, igcm_ho2)     = c(l,i_ho2)/dens(l) 
+         zycol(l, igcm_h2o2)    = c(l,i_h2o2)/dens(l)
+         zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
+         zycol(l, igcm_n)       = c(l,i_n)/dens(l)
+         zycol(l, igcm_n2d)     = c(l,i_n2d)/dens(l)
+         zycol(l, igcm_no)      = c(l,i_no)/dens(l)
+         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
+         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
+      end do 
+
+      end subroutine chimtogcm
+
+end subroutine photochemistry
+
Index: unk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90	(revision 2006)
+++ 	(revision )
@@ -1,1853 +1,0 @@
-!*****************************************************************
-!
-!     Photochemical routine 
-!
-!     Author: Franck Lefevre
-!     ------
-!
-!     Version: 27/04/2017
-!
-!     ASIS scheme : for details on the method see
-!     Cariolle et al., Geosci. Model Dev., 10, 1467-1485, 2017.
-!
-!*****************************************************************
-
-subroutine photochemistry_asis(nlayer, nq,                            &
-                          ig, lswitch, zycol, sza, ptimestep, press,  &
-                          temp, dens, zmmean, dist_sol, surfdust1d,   &
-                          surfice1d, jo3, tau, iter)
-
-implicit none
-
-#include "chimiedata.h" 
-#include "callkeys.h"
-
-!===================================================================
-!     inputs:
-!===================================================================
-
-integer, intent(in) :: nlayer ! number of atmospheric layers
-integer, intent(in) :: nq     ! number of tracers
-integer :: ig                 ! grid point index
-      
-real :: sza                   ! solar zenith angle (deg)
-real :: ptimestep             ! physics timestep (s)
-real :: press(nlayer)         ! pressure (hpa)
-real :: temp(nlayer)          ! temperature (k)
-real :: dens(nlayer)          ! density (cm-3)
-real :: zmmean(nlayer)        ! mean molar mass (g/mole)
-real :: dist_sol              ! sun distance (au) 
-real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
-real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
-real :: tau                   ! optical depth at 7 hpa
-
-!===================================================================
-!     input/output:
-!===================================================================
-      
-real :: zycol(nlayer,nq)       ! chemical species volume mixing ratio
-
-!===================================================================
-!     output:
-!===================================================================
-      
-integer :: iter(nlayer)        ! iteration counter
-real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
-
-!===================================================================
-!     local:
-!===================================================================
-
-integer :: phychemrat         ! (physical timestep)/(nominal chemical timestep)
-integer :: j_o3_o1d, ilev
-integer :: iesp, nesp
-integer :: lswitch
-
-logical, save :: firstcall = .true.
-
-parameter (nesp = 17)         ! number of species in the chemistry
-
-! tracer indexes in the chemistry:
-
-integer,parameter :: i_co2  =  1
-integer,parameter :: i_co   =  2
-integer,parameter :: i_o    =  3
-integer,parameter :: i_o1d  =  4
-integer,parameter :: i_o2   =  5
-integer,parameter :: i_o3   =  6
-integer,parameter :: i_h    =  7
-integer,parameter :: i_h2   =  8
-integer,parameter :: i_oh   =  9
-integer,parameter :: i_ho2  = 10
-integer,parameter :: i_h2o2 = 11
-integer,parameter :: i_h2o  = 12
-integer,parameter :: i_n    = 13
-integer,parameter :: i_n2d  = 14
-integer,parameter :: i_no   = 15
-integer,parameter :: i_no2  = 16
-integer,parameter :: i_n2   = 17
-
-real :: ctimestep           ! standard timestep for the chemistry (s) 
-real :: dt_guess            ! first-guess timestep (s) 
-real :: dt_corrected        ! corrected timestep (s) 
-real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
-real :: time                ! internal time (between 0 and ptimestep, in s)
-
-real, dimension(nlayer,nesp)            :: rm   ! mixing ratios 
-real (kind = 8), dimension(nesp)        :: cold ! number densities at previous timestep (molecule.cm-3) 
-real (kind = 8), dimension(nlayer,nesp) :: c    ! number densities at current timestep (molecule.cm-3) 
-real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3) 
- 
-! reaction rates
- 
-real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
-real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
-real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
-logical :: hetero_dust, hetero_ice
-
-! matrix
-
-real (kind = 8), dimension(nesp,nesp) :: mat, mat1
-integer, dimension(nesp)              :: indx
-integer                               :: code
-
-! production and loss terms (for first-guess solution only)
-
-real (kind = 8), dimension(nesp) :: prod, loss
-
-!===================================================================
-!     initialisation of the reaction indexes
-!===================================================================
-
-if (firstcall) then
-   print*,'photochemistry: initialize indexes'
-   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-               i_n, i_n2d, i_no, i_no2, i_n2)
-   firstcall = .false.
-end if
-
-!===================================================================
-!     initialisation of mixing ratios and densities       
-!===================================================================
-
-call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-               i_n, i_n2d, i_no, i_no2, i_n2, dens, rm, c) 
-
-!===================================================================
-!     interpolation of photolysis rates in the lookup table      
-!===================================================================
-
-call photolysis_asis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & 
-                     rm(:,i_co2), rm(:,i_o3), v_phot)
-
-! save o3 photolysis for output
-
-j_o3_o1d = 5
-jo3(:) = v_phot(:,j_o3_o1d)
-
-!===================================================================
-!     reaction rates                                     
-!===================================================================
-!     switches for heterogeneous chemistry
-!     hetero_ice  : reactions on ice clouds
-!     hetero_dust : reactions on dust    
-!===================================================================
-
-hetero_dust = .false.
-hetero_ice  = .true.
-
-call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2), &
-                   press, temp, hetero_dust, hetero_ice,         &
-                   surfdust1d, surfice1d, v_phot, v_3, v_4)
-
-!===================================================================
-!     ctimestep : standard chemical timestep (s), defined as 
-!                 the fraction phychemrat of the physical timestep                           
-!===================================================================
-
-phychemrat = 1
-
-ctimestep = ptimestep/real(phychemrat)
-
-!print*, "ptimestep  = ", ptimestep
-!print*, "phychemrat = ", phychemrat
-!print*, "ctimestep  = ", ctimestep
-!stop
-
-!===================================================================
-!     loop over levels         
-!===================================================================
-
-do ilev = 1,lswitch - 1
-
-!  initializations
-
-   time = 0.
-   iter(ilev) = 0
-   dt_guess = ctimestep
-   cold(:) = c(ilev,:)
-
-!  internal loop for the chemistry
-
-   do while (time < ptimestep)
-
-   iter(ilev) = iter(ilev) + 1    ! iteration counter
-  
-!  first-guess: fill matrix
-
-   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
-
-!  adaptative evaluation of the sub time step
-
-   call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:),  &
-                  mat1, prod, loss, dens(ilev))
-
-   if (time + dt_corrected > ptimestep) then
-      dt_corrected = ptimestep - time
-   end if
-
-!  if (dt_corrected /= dt_guess) then  ! the timestep has been modified
-
-!  form the matrix identity + mat*dt_corrected
-
-   mat(:,:) = mat1(:,:)*dt_corrected
-   do iesp = 1,nesp
-      mat(iesp,iesp) = 1. + mat(iesp,iesp)
-   end do
-
-!  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
-     
-   cnew(:) = c(ilev,:)
-
-#ifdef LAPACK
-   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
-#else
-   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-   stop
-#endif
-
-!  end if
-
-!  eliminate small values
-
-   where (cnew(:)/dens(ilev) < 1.e-30)
-      cnew(:) = 0.
-   end where
-
-!  update concentrations
-
-   cold(:)   = c(ilev,:)
-   c(ilev,:) = cnew(:)
-   cnew(:)   = 0.
-
-!  increment internal time
-
-   time = time + dt_corrected
-   dt_guess = dt_corrected     ! first-guess timestep for next iteration
-
-   end do ! while (time < ptimestep)
-
-end do ! ilev
-
-!===================================================================
-!     save chemical species for the gcm       
-!===================================================================
-
-call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
-               i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
-
-contains
-
-!================================================================
-
- subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
-                      prod, loss, dens)
-
-!================================================================
-! iterative evaluation of the appropriate time step dtnew
-! according to curvature criterion based on
-! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
-! with r = (tn - tn-1)/(tn+1 - tn)
-!================================================================
-
-implicit none
-
-! input
-
-integer :: nesp  ! number of species in the chemistry
-
-real :: dtold, ctimestep
-real (kind = 8), dimension(nesp)      :: cold, ccur
-real (kind = 8), dimension(nesp,nesp) :: mat1
-real (kind = 8), dimension(nesp)      :: prod, loss
-real                                  :: dens
-
-! output
-
-real :: dtnew
-
-! local
-
-real (kind = 8), dimension(nesp)      :: cnew
-real (kind = 8), dimension(nesp,nesp) :: mat
-real (kind = 8) :: atol, ratio, e, es, coef
-
-integer                  :: code, iesp, iter
-integer, dimension(nesp) :: indx
-
-real :: dttest
-
-! parameters
-
-real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
-real (kind = 8), parameter :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
-real (kind = 8), parameter :: rtol    = 1./0.05   ! 1/rtol recommended value : 0.1-0.02
-integer,         parameter :: niter   = 3        ! number of iterations
-real (kind = 8), parameter :: coefmax = 2.
-real (kind = 8), parameter :: coefmin = 0.1 
-logical                    :: fast_guess = .true.
-
-dttest = dtold   ! dttest = dtold = dt_guess
-
-atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
-
-do iter = 1,niter
-
-if (fast_guess) then
-
-! first guess : fast semi-implicit method
-
-   do iesp = 1, nesp
-      cnew(iesp) = (ccur(iesp) + prod(iesp)*dttest)/(1. + loss(iesp)*dttest)
-   end do
-
-else
-
-! first guess : form the matrix identity + mat*dt_guess
-
-   mat(:,:) = mat1(:,:)*dttest
-   do iesp = 1,nesp
-      mat(iesp,iesp) = 1. + mat(iesp,iesp)
-   end do
-
-! form right-hand side (RHS) of the system
-
-   cnew(:) = ccur(:)
-
-! solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
-
-#ifdef LAPACK
-      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
-#else
-   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-   stop
-#endif
-
-end if
-
-! ratio old/new subtimestep
-
-ratio = dtold/dttest
-
-! e : local error indicatocitr
-
-e = 0.
-
-do iesp = 1,nesp
-   es = 2.*abs((ratio*cnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp))   &
-         /(1. + ratio)/max(ccur(iesp),atol))
-
-   if (es > e) then
-      e = es
-   end if
-end do
-e = rtol*e
-
-! timestep correction
-
-coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
-
-dttest = max(dtmin,dttest*coef)
-dttest = min(ctimestep,dttest)
-
-end do ! iter
-
-! new timestep
-
-dtnew = dttest
-
-end subroutine define_dt
-
-!======================================================================
-
- subroutine reactionrates(nlayer,                               &
-                          lswitch, dens, co2, o2, press, t,     &
-                          hetero_dust, hetero_ice,              &
-                          surfdust1d, surfice1d,                &
-                          v_phot, v_3, v_4)
- 
-!================================================================
-! compute reaction rates                                        !
-!----------------------------------------------------------------
-! reaction               type                array              !
-!----------------------------------------------------------------
-! A + B    --> C + D     bimolecular         v_4                !
-! A + A    --> B + C     quadratic           v_3                !
-! A + C    --> B + C     quenching           v_phot             !
-! A + ice  --> B + C     heterogeneous       v_phot             !
-!================================================================
-
-use comcstfi_h
-
-implicit none
-
-#include "chimiedata.h"
-
-!----------------------------------------------------------------------
-!     input
-!----------------------------------------------------------------------
-
-integer, intent(in)     :: nlayer            ! number of atmospheric layers
-integer                 :: lswitch           ! interface level between lower
-                                             ! atmosphere and thermosphere chemistries
-real, dimension(nlayer) :: dens              ! total number density (molecule.cm-3)
-real, dimension(nlayer) :: press             ! pressure (hPa)
-real, dimension(nlayer) :: t                 ! temperature (K)
-real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
-real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
-real (kind = 8), dimension(nlayer) :: co2    ! co2 number density (molecule.cm-3)
-real (kind = 8), dimension(nlayer) :: o2     ! o2 number density (molecule.cm-3)
-logical :: hetero_dust, hetero_ice           ! switches for heterogeneous chemistry
-
-!----------------------------------------------------------------------
-!     output
-!----------------------------------------------------------------------
-
-real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
-real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
-real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
-
-!----------------------------------------------------------------------
-!     local
-!----------------------------------------------------------------------
-
-integer          :: ilev
-integer          :: nb_phot, nb_reaction_3, nb_reaction_4
-real :: ak0, ak1, xpo, rate
-real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
-real, dimension(nlayer) :: deq
-real, dimension(nlayer) :: a001, a002, a003,                           &
-                             b001, b002, b003, b004, b005, b006, b007,   &
-                             b008, b009,                                 &
-                             c001, c002, c003, c004, c005, c006, c007,   &
-                             c008, c009, c010, c011, c012, c013, c014,   &
-                             c015, c016, c017, c018,                     &
-                             d001, d002, d003, d004, d005, d006, d007,   &
-                             d008, d009,                                 &
-                             e001, e002,                                 &
-                             h001, h002, h003, h004, h005
-
-!----------------------------------------------------------------------
-!     initialisation
-!----------------------------------------------------------------------
-
-      nb_phot       = 11       ! jmars.20140930 reduit de 13 a 11
-      nb_reaction_3 = 0
-      nb_reaction_4 = 0
-
-!----------------------------------------------------------------------
-!     oxygen reactions 
-!----------------------------------------------------------------------
-
-!---  a001: o + o2 + co2 -> o3 + co2
-
-!     jpl 2003
-!
-!     co2/n2 efficiency as a third body = 2.075
-!     from sehested et al., j. geophys. res., 100, 1995.
-
-      a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = a001(:)
-
-!---  a002: o + o + co2 -> o2 + co2
-
-!     Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
-
-!     a002(:) = 2.5*5.2e-35*exp(900./t(:))*dens(:)
-
-!     Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
-
-!     a002(:) = 1.2e-32*(300./t(:))**(2.0)*dens(:)  ! yung expression
-      a002(:) = 2.5*9.46e-34*exp(485./t(:))*dens(:) ! nist expression
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = a002(:)
-
-!---  a003: o + o3 -> o2 + o2
-
-!     jpl 2003
-
-      a003(:) = 8.0e-12*exp(-2060./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = a003(:)
-
-!----------------------------------------------------------------------
-!     o(1d) reactions
-!----------------------------------------------------------------------
-
-!---  b001: o(1d) + co2  -> o + co2
-
-!     jpl 2006
-
-      b001(:) = 7.5e-11*exp(115./t(:))
-   
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = b001(:)*co2(:)
-
-!---  b002: o(1d) + h2o  -> oh + oh
-
-!     jpl 2006
- 
-      b002(:) = 1.63e-10*exp(60./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b002(:)
-
-!---  b003: o(1d) + h2  -> oh + h
-
-!     jpl 2011
-
-      b003(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b003(:)
-
-!---  b004: o(1d) + o2  -> o + o2
-
-!     jpl 2006
-
-      b004(:) = 3.3e-11*exp(55./t(:))
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = b004(:)*o2(:)
-    
-!---  b005: o(1d) + o3  -> o2 + o2
-
-!     jpl 2003
-
-      b005(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b005(:)
-    
-!---  b006: o(1d) + o3  -> o2 + o + o
-
-!     jpl 2003
-
-      b006(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b006(:)
-    
-!---  b007: o(1d) + ch4 -> ch3 + oh
-
-!     jpl 2003
-
-      b007(:) = 1.5e-10*0.75
-
-!---  b008: o(1d) + ch4 -> ch3o + h
-
-!     jpl 2003
-
-      b008(:) = 1.5e-10*0.20
-!
-!---  b009: o(1d) + ch4 -> ch2o + h2
-
-!     jpl 2003
-
-      b009(:) = 1.5e-10*0.05
-
-!----------------------------------------------------------------------
-!     hydrogen reactions
-!----------------------------------------------------------------------
-
-!---  c001: o + ho2 -> oh + o2
-
-!     jpl 2003
-
-      c001(:) = 3.0e-11*exp(200./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c001(:)
-
-!---  c002: o + oh -> o2 + h
-
-!     jpl 2011
-
-      c002(:) = 1.8e-11*exp(180./t(:))
-
-!     robertson and smith, j. chem. phys. a 110, 6673, 2006
-
-!     c002(:) = 11.2e-11*t(:)**(-0.32)*exp(177./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c002(:)
-
-!---  c003: h + o3 -> oh + o2
-
-!     jpl 2003
-
-      c003(:) = 1.4e-10*exp(-470./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c003(:)
-
-!---  c004: h + ho2 -> oh + oh
-
-!     jpl 2006
-
-      c004(:) = 7.2e-11
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c004(:)
-
-!---  c005: h + ho2 -> h2 + o2
-
-!     jpl 2006
-
-      c005(:) = 6.9e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c005(:)
-
-!---  c006: h + ho2 -> h2o + o
-
-!     jpl 2006
-
-      c006(:) = 1.6e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c006(:)
-
-!---  c007: oh + ho2 -> h2o + o2
-
-!     jpl 2003
-
-!     canty et al., grl, 2006 suggest to increase this rate
-!     by 20%. not done here.
-
-      c007(:) = 4.8e-11*exp(250./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c007(:)
-
-!---  c008: ho2 + ho2 -> h2o2 + o2
-
-!     jpl 2006
-
-!     c008(:) = 3.5e-13*exp(430./t(:))
-
-!     christensen et al., grl, 13, 2002
-
-      c008(:) = 1.5e-12*exp(19./t(:))
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c008(:)
-
-!---  c009: oh + h2o2 -> h2o + ho2
-
-!     jpl 2006
-
-      c009(:) = 1.8e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c009(:)
-
-!---  c010: oh + h2 -> h2o + h
-
-!     jpl 2006
-
-      c010(:) = 2.8e-12*exp(-1800./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c010(:)
-
-!---  c011: h + o2 + co2 -> ho2 + co2
-
-!     jpl 2011
-!     co2/n2 efficiency as a third body = 2.4
-!     from ashman and haynes, 27th symposium on combustion, 1998.
-
-      do ilev = 1,lswitch-1
-         ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3)
-         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
-
-         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
-         c011(ilev) = rate*0.6**xpo
-      end do
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c011(:)
-
-!---  c012: o + h2o2 -> oh + ho2
-
-!     jpl 2003
-
-      c012(:) = 1.4e-12*exp(-2000./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c012(:)
-
-!---  c013: oh + oh -> h2o + o
-
-!     jpl 2006
-
-      c013(:) = 1.8e-12
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c013(:)
-
-!---  c014: oh + o3 -> ho2 + o2
-
-!     jpl 2003
-
-      c014(:) = 1.7e-12*exp(-940./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c014(:)
-
-!---  c015: ho2 + o3 -> oh + o2 + o2
-
-!     jpl 2003
-
-      c015(:) = 1.0e-14*exp(-490./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c015(:)
-
-!---  c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
-
-!     jpl 2011
-
-      c016(:) = 2.5*2.1e-33*exp(920./t(:))*dens(:)
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c016(:)
-
-!---  c017: oh + oh + co2 -> h2o2 + co2
-
-!     jpl 2003
-
-      do ilev = 1,lswitch-1
-         ak0 = 2.5*6.9e-31*(t(ilev)/300.)**(-1.0)
-         ak1 = 2.6e-11*(t(ilev)/300.)**(0.0)
-
-         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
-         c017(ilev) = rate*0.6**xpo
-      end do
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c017(:)
-
-!---  c018: h + h + co2 -> h2 + co2
-
-!     baulch et al., 2005
-
-      c018(:) = 2.5*1.8e-30*(t(:)**(-1.0))*dens(:)
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c018(:)
-
-!----------------------------------------------------------------------
-!     nitrogen reactions
-!----------------------------------------------------------------------
-
-!---  d001: no2 + o -> no + o2
-
-!     jpl 2006
-
-      d001(:) = 5.1e-12*exp(210./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d001(:)
-
-!---  d002: no + o3 -> no2 + o2
-
-!     jpl 2006
-
-      d002(:) = 3.0e-12*exp(-1500./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d002(:)
-
-!---  d003: no + ho2 -> no2 + oh
-
-!     jpl 2011
-
-      d003(:) = 3.3e-12*exp(270./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d003(:)
-
-!---  d004: n + no -> n2 + o
-
-!     jpl 2011
-
-      d004(:) = 2.1e-11*exp(100./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d004(:)
-
-!---  d005: n + o2 -> no + o
-
-!     jpl 2011
-
-      d005(:) = 1.5e-11*exp(-3600./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d005(:)
-
-!---  d006: no2 + h -> no + oh
-
-!     jpl 2011
-
-      d006(:) = 4.0e-10*exp(-340./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d006(:)
-
-!---  d007: n + o -> no
- 
-      d007(:) = 2.8e-17*(300./t(:))**0.5
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d007(:)
-
-!---  d008: n + ho2 -> no + oh
-
-!     brune et al., j. chem. phys., 87, 1983 
-
-      d008(:) = 2.19e-11
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d008(:)
-
-!---  d009: n + oh -> no + h
-
-!     atkinson et al., j. phys. chem. ref. data, 18, 881, 1989 
-
-      d009(:) = 3.8e-11*exp(85./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d009(:)
-
-!----------------------------------------------------------------------
-!     carbon reactions
-!----------------------------------------------------------------------
-
-!---  e001: oh + co -> co2 + h
-
-!     jpl 2003
-
-!     e001(:) = 1.5e-13*(1 + 0.6*press(:)/1013.)
-
-!     mccabe et al., grl, 28, 3135, 2001
-
-!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
-
-!     jpl 2006
-
-!     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
-!     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
-!     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
-!     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
-
-!     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
-!     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
-!     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
-!     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
-
-!     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
-
-!     joshi et al., 2006
-
-      do ilev = 1,lswitch-1
-         k1a0 = 1.34*2.5*dens(ilev)                                  &
-               *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
-               + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
-         k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
-              + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
-         k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
-                + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
-         x = k1a0/(k1ainf - k1b0)
-         y = k1b0/(k1ainf - k1b0)
-         fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
-            + exp(-t(ilev)/255.)
-         fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
-         k1a = k1a0*((1. + y)/(1. + x))*fx
-         k1b = k1b0*(1./(1.+x))*fx
-
-         e001(ilev) = k1a + k1b
-      end do
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = e001(:)
-
-!---  e002: o + co + m -> co2 + m
-
-!     tsang and hampson, 1986.
-
-      e002(:) = 2.5*6.5e-33*exp(-2184./t(:))*dens(:)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = e002(:)
-
-!----------------------------------------------------------------------
-!     heterogeneous chemistry 
-!----------------------------------------------------------------------
-
-      if (hetero_ice) then
-
-!        k = (surface*v*gamma)/4 (s-1)
-!        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
- 
-!---     h001: ho2 + ice -> products
- 
-!        cooper and abbatt, 1996: gamma = 0.025
-      
-         gam = 0.025
-         h001(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
- 
-!        h002: oh + ice -> products
- 
-!        cooper and abbatt, 1996: gamma = 0.03
- 
-         gam = 0.03
-         h002(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(17.e-3*pi))*gam/4.
-
-!---     h003: h2o2 + ice -> products
- 
-!        gamma = 0.    test value
- 
-         gam = 0.
-         h003(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
-      else
-         h001(:) = 0.
-         h002(:) = 0.
-         h003(:) = 0.
-      end if
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h001(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h002(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h003(:)
-
-      if (hetero_dust) then
- 
-!---     h004: ho2 + dust -> products
- 
-!        jacob, 2000: gamma = 0.2
-!        see dereus et al., atm. chem. phys., 2005
- 
-         gam = 0.2
-         h004(:) = surfdust1d(:)  &
-                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
- 
-!---     h005: h2o2 + dust -> products
- 
-!        gamma = 5.e-4
-!        see dereus et al., atm. chem. phys., 2005
- 
-         gam = 5.e-4
-         h005(:) = surfdust1d(:)  &
-                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
-      else
-         h004(:) = 0.
-         h005(:) = 0.
-      end if
- 
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h004(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h005(:)
-
-end subroutine reactionrates
-
-!======================================================================
-
- subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
-
-!======================================================================
-! filling of the jacobian matrix
-!======================================================================
-
-use types_asis
-
-implicit none
-
-#include "chimiedata.h"
-
-! input
-
-integer             :: ilev    ! level index
-integer             :: nesp    ! number of species in the chemistry
-integer, intent(in) :: nlayer  ! number of atmospheric layers
-
-real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
-real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
-real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
-real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
-
-! output
-
-real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
-real (kind = 8), dimension(nesp), intent(out)      :: prod, loss
-
-! local
-
-integer :: iesp
-integer :: ind_phot_2,ind_phot_4,ind_phot_6
-integer :: ind_3_2,ind_3_4,ind_3_6
-integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
-integer :: iphot,i3,i4
-
-real(kind = 8) :: eps, eps_4  ! implicit/explicit coefficient
-
-! initialisations 
-
-mat(:,:) = 0.
-prod(:)  = 0.
-loss(:)  = 0.
-
-! photodissociations
-! or reactions a + c -> b + c
-! or reactions a + ice -> b + c
-
-do iphot = 1,nb_phot_max
-
-  ind_phot_2 = indice_phot(iphot)%z2
-  ind_phot_4 = indice_phot(iphot)%z4
-  ind_phot_6 = indice_phot(iphot)%z6
-
-  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
-  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
-  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
-
-  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
-  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
-  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
-
-end do
-
-! reactions a + a -> b + c 
-
-do i3 = 1,nb_reaction_3_max
-
-  ind_3_2 = indice_3(i3)%z2
-  ind_3_4 = indice_3(i3)%z4
-  ind_3_6 = indice_3(i3)%z6
-
-  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
-  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)
-  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)
-
-  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
-  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
-  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
-
-end do
-
-! reactions a + b -> c + d 
-
-eps = 1.d-10
-
-do i4 = 1,nb_reaction_4_max
-
-  ind_4_2 = indice_4(i4)%z2
-  ind_4_4 = indice_4(i4)%z4
-  ind_4_6 = indice_4(i4)%z6
-  ind_4_8 = indice_4(i4)%z8
-
-  eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
-  eps_4 = min(eps_4,1.0)
-
-  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4) 
-  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
-  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
-  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)   
-  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
-  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
-  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
-  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)
-  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
-  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
-
-end do
-
-end subroutine fill_matrix
-
-!================================================================
-
- 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)
-
-!================================================================
-! set the "indice" arrays used to fill the jacobian matrix      !
-!----------------------------------------------------------------
-! reaction                                   type               !
-!----------------------------------------------------------------
-! A + hv   --> B + C     photolysis          indice_phot        ! 
-! A + B    --> C + D     bimolecular         indice_4           !
-! A + A    --> B + C     quadratic           indice_3           !
-! A + C    --> B + C     quenching           indice_phot        !
-! A + ice  --> B + C     heterogeneous       indice_phot        !
-!================================================================
-
-use types_asis
-
-implicit none
-
-#include "chimiedata.h"
-
-! input
-
-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
-
-! local
-
-integer :: nb_phot, nb_reaction_3, nb_reaction_4
-integer :: i_dummy
-
-allocate (indice_phot(nb_phot_max))
-allocate (indice_3(nb_reaction_3_max))
-allocate (indice_4(nb_reaction_4_max))
-
-i_dummy = 1
-
-nb_phot       = 0
-nb_reaction_3 = 0
-nb_reaction_4 = 0
-
-!===========================================================
-!      O2 + hv -> O + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      O2 + hv -> O + O(1D)
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d) 
-
-!===========================================================
-!      CO2 + hv -> CO + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o) 
-
-!===========================================================
-!      CO2 + hv -> CO + O(1D)
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d) 
-
-!===========================================================
-!      O3 + hv -> O2 + O(1D)
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d) 
-
-!===========================================================
-!      O3 + hv -> O2 + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o) 
-
-!===========================================================
-!      H2O + hv -> H + OH
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh) 
-
-!===========================================================
-!      H2O2 + hv -> OH + OH
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      HO2 + hv -> OH + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o) 
-
-!===========================================================
-!      NO + hv -> N + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o) 
-
-!===========================================================
-!      NO2 + hv -> NO + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o) 
-
-!===========================================================
-!      a001 : O + O2 + CO2 -> O3 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy) 
-
-!===========================================================
-!      a002 : O + O + CO2 -> O2 + CO2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      a003 : O + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b001 : O(1D) + CO2 -> O + CO2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b002 : O(1D) + H2O -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      b003 : O(1D) + H2 -> OH + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h) 
-
-!===========================================================
-!      b004 : O(1D) + O2 -> O + O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b005 : O(1D) + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b006 : O(1D) + O3 -> O2 + O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 
-
-!===========================================================
-!      c001 : O + HO2 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c002 : O + OH -> O2 + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h) 
-
-!===========================================================
-!      c003 : H + O3 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c004 : H + HO2 -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      c005 : H + HO2 -> H2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2) 
-
-!===========================================================
-!      c006 : H + HO2 -> H2O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c007 : OH + HO2 -> H2O + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2) 
-
-!===========================================================
-!      c008 : HO2 + HO2 -> H2O2 + O2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c009 : OH + H2O2 -> H2O + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2) 
-
-!===========================================================
-!      c010 : OH + H2 -> H2O + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h) 
-
-!===========================================================
-!      c011 : H + O2 + CO2 -> HO2 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy) 
-
-!===========================================================
-!      c012 : O + H2O2 -> OH + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2) 
-
-!===========================================================
-!      c013 : OH + OH -> H2O + O
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c014 : OH + O3 -> HO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2) 
-
-!===========================================================
-!      c015 : HO2 + O3 -> OH + O2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2) 
-
-!===========================================================
-!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c017 : OH + OH + CO2 -> H2O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy) 
-
-!===========================================================
-!      c018 : H + H + CO2 -> H2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 
-
-!===========================================================
-!      d001 : NO2 + O -> NO + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2) 
-
-!===========================================================
-!      d002 : NO + O3 -> NO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2) 
-
-!===========================================================
-!      d003 : NO + HO2 -> NO2 + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh) 
-
-!===========================================================
-!      d004 : N + NO -> N2 + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o) 
-
-!===========================================================
-!      d005 : N + O2 -> NO + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o) 
-
-!===========================================================
-!      d006 : NO2 + H -> NO + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d007 : N + O -> NO
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy) 
-
-!===========================================================
-!      d008 : N + HO2 -> NO + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d009 : N + OH -> NO + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h) 
-
-!===========================================================
-!      e001 : CO + OH -> CO2 + H 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h) 
-
-!===========================================================
-!      e002 : CO + O + M -> CO2 + M 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 
-
-!===========================================================
-!      h001: HO2 + ice -> products
-!            treated as
-!            HO2 -> 0.5 H2O + 0.75 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
-
-!===========================================================
-!      h002: OH + ice -> products
-!            treated as
-!            OH -> 0.5 H2O + 0.25 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2) 
-
-!===========================================================
-!      h003: H2O2 + ice -> products
-!            treated as
-!            H2O2 -> H2O + 0.5 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
-
-!===========================================================
-!      h004: HO2 + dust -> products
-!            treated as
-!            HO2 -> 0.5 H2O + 0.75 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
-
-!===========================================================
-!      h005: H2O2 + dust -> products
-!            treated as
-!            H2O2 -> H2O + 0.5 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
-
-!===========================================================
-!  check dimensions 
-!===========================================================
-
-print*, 'nb_phot       = ', nb_phot
-print*, 'nb_reaction_4 = ', nb_reaction_4
-print*, 'nb_reaction_3 = ', nb_reaction_3
-
-if ((nb_phot /= nb_phot_max)             .or.  &
-    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
-    (nb_reaction_4 /= nb_reaction_4_max)) then
-   print*, 'wrong dimensions in indice' 
-   stop
-end if  
-
-end subroutine indice
-
-!*****************************************************************
-
-      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
-                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
-                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
-                           i_n, i_n2d, i_no, i_no2, i_n2,            &
-                           dens, rm, c) 
-
-!*****************************************************************
-
-      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
-     &                      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
-
-      implicit none
-
-#include "callkeys.h"
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     input:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer, intent(in) :: nq     ! number of tracers in the gcm
-      integer :: nesp               ! number of species in the chemistry
-      integer :: lswitch            ! interface level between chemistries
-
-      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
-
-      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
-      real :: dens(nlayer)          ! total number density (molecule.cm-3) 
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     output:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real, dimension(nlayer,nesp)            :: rm ! volume mixing ratios 
-      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 
-      
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     local:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer      :: l, iesp
-      logical,save :: firstcall = .true.
-      
-      
-!     first call initializations
-
-      if (firstcall) then
-
-!       identify the indexes of the tracers we need
-
-         if (igcm_co2 == 0) then
-            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
-            stop
-         endif
-         if (igcm_co == 0) then
-            write(*,*) "gcmtochim: Error; no CO tracer !!!"
-            stop
-         end if
-         if (igcm_o == 0) then
-            write(*,*) "gcmtochim: Error; no O tracer !!!"
-            stop
-         end if
-         if (igcm_o1d == 0) then
-            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
-            stop
-         end if
-         if (igcm_o2 == 0) then
-            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
-            stop
-         end if
-         if (igcm_o3 == 0) then
-            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
-            stop
-         end if
-         if (igcm_h == 0) then
-            write(*,*) "gcmtochim: Error; no H tracer !!!"
-            stop
-         end if
-         if (igcm_h2 == 0) then
-            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
-            stop
-         end if
-         if (igcm_oh == 0) then
-            write(*,*) "gcmtochim: Error; no OH tracer !!!"
-            stop
-         end if
-         if (igcm_ho2 == 0) then
-            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
-            stop
-         end if
-         if (igcm_h2o2 == 0) then
-            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
-            stop
-         end if
-         if (igcm_n == 0) then
-            write(*,*) "gcmtochim: Error; no N tracer !!!"
-            stop
-         end if
-         if (igcm_n2d == 0) then
-            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
-            stop
-         end if
-         if (igcm_no == 0) then
-            write(*,*) "gcmtochim: Error; no NO tracer !!!"
-            stop
-         end if
-         if (igcm_no2 == 0) then
-            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
-            stop
-         end if
-         if (igcm_n2 == 0) then
-            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
-            stop
-         end if
-         if (igcm_h2o_vap == 0) then
-            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
-            stop
-         end if
-         firstcall = .false.
-      end if ! of if (firstcall)
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     initialise mixing ratios
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         rm(l,i_co2)  = zycol(l, igcm_co2)
-         rm(l,i_co)   = zycol(l, igcm_co)
-         rm(l,i_o)    = zycol(l, igcm_o)
-         rm(l,i_o1d)  = zycol(l, igcm_o1d)
-         rm(l,i_o2)   = zycol(l, igcm_o2)
-         rm(l,i_o3)   = zycol(l, igcm_o3)
-         rm(l,i_h)    = zycol(l, igcm_h)
-         rm(l,i_h2)   = zycol(l, igcm_h2)
-         rm(l,i_oh)   = zycol(l, igcm_oh)
-         rm(l,i_ho2)  = zycol(l, igcm_ho2)
-         rm(l,i_h2o2) = zycol(l, igcm_h2o2)
-         rm(l,i_h2o)  = zycol(l, igcm_h2o_vap)
-         rm(l,i_n)    = zycol(l, igcm_n)
-         rm(l,i_n2d)  = zycol(l, igcm_n2d)
-         rm(l,i_no)   = zycol(l, igcm_no)
-         rm(l,i_no2)  = zycol(l, igcm_no2)
-         rm(l,i_n2)   = zycol(l, igcm_n2)
-      end do 
-
-      where (rm(:,:) < 1.e-30)
-         rm(:,:) = 0.
-      end where
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     initialise number densities
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do iesp = 1,nesp
-         do l = 1,lswitch-1
-            c(l,iesp) = rm(l,iesp)*dens(l)
-         end do
-      end do
-
-      end subroutine gcmtochim
-
-!*****************************************************************
- 
-      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,         &
-                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
-                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
-                           i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
- 
-!*****************************************************************
- 
-      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,          &
-                            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
-
-      implicit none
-
-#include "callkeys.h"
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     inputs:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
- 
-      integer, intent(in) :: nlayer  ! number of atmospheric layers
-      integer, intent(in) :: nq      ! number of tracers in the gcm
-      integer :: nesp                ! number of species in the chemistry
-      integer :: lswitch             ! interface level between chemistries
-      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
-
-      real :: dens(nlayer)     ! total number density (molecule.cm-3) 
-      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     output:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-       
-      real zycol(nlayer,nq)  ! volume mixing ratios in the gcm
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     local:
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer l
-      
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!     save mixing ratios for the gcm
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-         zycol(l, igcm_co2)     = c(l,i_co2)/dens(l) 
-         zycol(l, igcm_co)      = c(l,i_co)/dens(l) 
-         zycol(l, igcm_o)       = c(l,i_o)/dens(l) 
-         zycol(l, igcm_o1d)     = c(l,i_o1d)/dens(l)
-         zycol(l, igcm_o2)      = c(l,i_o2)/dens(l) 
-         zycol(l, igcm_o3)      = c(l,i_o3)/dens(l) 
-         zycol(l, igcm_h)       = c(l,i_h)/dens(l)  
-         zycol(l, igcm_h2)      = c(l,i_h2)/dens(l) 
-         zycol(l, igcm_oh)      = c(l,i_oh)/dens(l) 
-         zycol(l, igcm_ho2)     = c(l,i_ho2)/dens(l) 
-         zycol(l, igcm_h2o2)    = c(l,i_h2o2)/dens(l)
-         zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
-         zycol(l, igcm_n)       = c(l,i_n)/dens(l)
-         zycol(l, igcm_n2d)     = c(l,i_n2d)/dens(l)
-         zycol(l, igcm_no)      = c(l,i_no)/dens(l)
-         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
-         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
-      end do 
-
-      end subroutine chimtogcm
-
-end subroutine photochemistry_asis
-
Index: /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 2006)
+++ /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 2007)
@@ -3,5 +3,5 @@
       subroutine photolysis(nlayer,                              &
                             lswitch, press, temp, sza, tauref,   &
-                            dist_sol, rmco2, rmo3, j)
+                            zmmean, dist_sol, rmco2, rmo3, v_phot)
 
 !==========================================================================
@@ -23,4 +23,5 @@
       real :: sza                   ! solar zenith angle (deg)
       real :: tauref                ! optical depth at 7 hpa
+      real :: zmmean(nlayer)        ! mean molecular mass (g)
       real :: dist_sol              ! sun distance (AU) 
       real :: rmco2(nlayer)         ! co2 mixing ratio
@@ -31,5 +32,5 @@
 !==========================================================================
 
-      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
+      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot
 
 !==========================================================================
@@ -39,4 +40,12 @@
       integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
                  iozo, isza, itau, it, l
+
+      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d,        &
+                 j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2,      &
+                 j_hno3, j_hno4,                                        &
+                 j_ch4_ch3_h, j_ch4_1ch2_h2, j_ch4_3ch2_h_h,            &
+                 j_ch4_ch_h2_h, j_ch3o2h, j_ch2o_hco, j_ch2o_co,        &
+                 j_ch3oh, j_c2h6, j_hcl, j_hocl, j_clo, j_so2, j_so,    &
+                 j_h2s, j_so3
 
       real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
@@ -49,4 +58,5 @@
       real :: ratio_o3(nlayer)
       real :: tau
+      real :: j(nlayer,nd)
 
 !==========================================================================
@@ -105,5 +115,5 @@
 
       col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
-                       /(mugaz*g)
+                       /(zmmean(lswitch-1)*g)
 
 !     ozone column at model top
@@ -116,8 +126,8 @@
          dp = (press(l) - press(l+1))*100.
          col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
-                             *6.022e22*dp/(mugaz*g)
+                             *6.022e22*dp/(zmmean(l)*g)
          col(l) = min(col(l), colairtab(1))
          colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
-                                 *6.022e22*dp/(mugaz*g)
+                                 *6.022e22*dp/(zmmean(l)*g)
       end do
 
@@ -295,4 +305,74 @@
       end if
 
+! photodissociation rates numbering in the lookup table
+
+! jmars.20140930
+
+      j_o2_o         =  1      ! o2 + hv     -> o + o
+      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
+      j_co2_o        =  3      ! co2 + hv    -> co + o
+      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
+      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
+      j_o3_o         =  6      ! o3 + hv     -> o2 + o
+      j_h2o          =  7      ! h2o + hv    -> h + oh
+      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
+      j_ho2          =  9      ! ho2 + hv    -> oh + o
+      j_no           =  10     ! no + hv     -> n + o
+      j_no2          =  11     ! no2 + hv    -> no + o
+      j_hno3         =  12     ! hno3 + hv   -> no2 + oh
+      j_hno4         =  13     ! hno4 + hv   -> no2 + ho2
+
+! jmars.20111014
+
+!     j_o2_o         =  1      ! o2 + hv     -> o + o
+!     j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
+!     j_co2_o        =  3      ! co2 + hv    -> co + o
+!     j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
+!     j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
+!     j_o3_o         =  6      ! o3 + hv     -> o2 + o
+!     j_h2o          =  7      ! h2o + hv    -> h + oh
+!     j_hdo          =  8      ! hdo + hv    -> d + oh
+!     j_h2o2         =  9      ! h2o2 + hv   -> oh + oh
+!     j_ho2          =  10     ! ho2 + hv    -> oh + o
+!     j_no2          =  11     ! no2 + hv    -> no + o
+!     j_ch4_ch3_h    =  12     ! ch4 + hv    -> ch3 + h
+!     j_ch4_1ch2_h2  =  13     ! ch4 + hv    -> 1ch2 + h2
+!     j_ch4_3ch2_h_h =  14     ! ch4 + hv    -> 3ch2 + h + h
+!     j_ch4_ch_h2_h  =  15     ! ch4 + hv    -> ch + h2 + h
+!     j_ch3o2h       =  16     ! ch3o2h + hv -> ch3o + oh
+!     j_ch2o_hco     =  17     ! ch2o + hv   -> h + hco
+!     j_ch2o_co      =  18     ! ch2o + hv   -> h2 + co
+!     j_ch3oh        =  19     ! ch3oh + hv  -> ch3o + h
+!     j_c2h6         =  20     ! c2h6 + hv   -> products
+!     j_hcl          =  21     ! hcl + hv    -> h + cl
+!     j_hocl         =  22     ! hocl + hv   -> oh + cl
+!     j_clo          =  23     ! clo + hv    -> o + cl
+!     j_so2          =  24     ! so2 + hv    -> so + o
+!     j_so           =  25     ! so + hv     -> s + o
+!     j_h2s          =  26     ! h2s + hv    -> hs + s
+!     j_so3          =  27     ! so2 + hv    -> so2 + o
+!     j_hno3         =  28     ! hno3 + hv   -> oh + no2
+!     j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
+
+! fill v_phot array
+
+      v_phot(:,:) = 0.
+
+      do l = 1,lswitch-1
+         v_phot(l, 1) = j(l,j_o2_o)
+         v_phot(l, 2) = j(l,j_o2_o1d)
+         v_phot(l, 3) = j(l,j_co2_o)
+         v_phot(l, 4) = j(l,j_co2_o1d)
+         v_phot(l, 5) = j(l,j_o3_o1d)
+         v_phot(l, 6) = j(l,j_o3_o)
+         v_phot(l, 7) = j(l,j_h2o)
+         v_phot(l, 8) = j(l,j_h2o2)
+         v_phot(l, 9) = j(l,j_ho2)
+         v_phot(l,10) = j(l,j_no)
+         v_phot(l,11) = j(l,j_no2)
+      end do
+
       return
       end
+
+!*****************************************************************
Index: unk/LMDZ.MARS/libf/aeronomars/photolysis_asis.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photolysis_asis.F90	(revision 2006)
+++ 	(revision )
@@ -1,378 +1,0 @@
-!==========================================================================
-
-      subroutine photolysis_asis(nlayer,                              &
-                                 lswitch, press, temp, sza, tauref,   &
-                                 zmmean, dist_sol, rmco2, rmo3, v_phot)
-
-!==========================================================================
-
-      use comcstfi_h
-
-      implicit none
-
-#include "chimiedata.h" 
-
-!==========================================================================
-!     input:
-!==========================================================================
-       
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer :: lswitch            ! interface level between chemistries
-      real :: press(nlayer)         ! pressure (hPa)
-      real :: temp(nlayer)          ! temperature (K)
-      real :: sza                   ! solar zenith angle (deg)
-      real :: tauref                ! optical depth at 7 hpa
-      real :: zmmean(nlayer)        ! mean molecular mass (g)
-      real :: dist_sol              ! sun distance (AU) 
-      real :: rmco2(nlayer)         ! co2 mixing ratio
-      real :: rmo3(nlayer)          ! ozone mixing ratio
-
-!==========================================================================
-!     output: interpolated photodissociation rates (s-1)
-!==========================================================================
-
-      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot
-
-!==========================================================================
-!     local:
-!==========================================================================
-
-      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
-                 iozo, isza, itau, it, l
-
-      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d,        &
-                 j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2,      &
-                 j_hno3, j_hno4,                                        &
-                 j_ch4_ch3_h, j_ch4_1ch2_h2, j_ch4_3ch2_h_h,            &
-                 j_ch4_ch_h2_h, j_ch3o2h, j_ch2o_hco, j_ch2o_co,        &
-                 j_ch3oh, j_c2h6, j_hcl, j_hocl, j_clo, j_so2, j_so,    &
-                 j_h2s, j_so3
-
-      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
-      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
-      real :: poids(2,2,2,2,2)            ! 5D interpolation weights 
-      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
-      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
-      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau
-      real :: colo3min, dp, coef
-      real :: ratio_o3(nlayer)
-      real :: tau
-      real :: j(nlayer,nd)
-
-!==========================================================================
-!     day/night criterion
-!==========================================================================
-
-      if (sza <= 95.) then 
-
-!==========================================================================
-!     temperatures at 1.9 hPa in lookup table
-!==========================================================================
-      
-      table_temp(1) = 226.2
-      table_temp(2) = 206.2
-      table_temp(3) = 186.2
-      table_temp(4) = 169.8
-
-!==========================================================================
-!     interpolation in solar zenith angle
-!==========================================================================
- 
-      indsza = nsza - 1
-      do isza = 1,nsza
-         if (szatab(isza) >= sza) then
-            indsza = isza - 1
-            indsza = max(indsza, 1)
-            exit
-         end if
-      end do
-      cisza = (sza - szatab(indsza))  &
-             /(szatab(indsza + 1) - szatab(indsza))
-
-!==========================================================================
-!     interpolation in dust (tau)
-!==========================================================================
-
-      tau = min(tauref, tautab(ntau))
-      tau = max(tau, tautab(1))
-
-      indtau = ntau - 1
-      do itau = 1,ntau
-         if (tautab(itau) >= tau) then
-            indtau = itau - 1
-            indtau = max(indtau, 1)
-            exit
-         end if
-      end do
-      citau = (tau - tautab(indtau))     &
-             /(tautab(indtau + 1) - tautab(indtau))
-
-!==========================================================================
-!     co2 and ozone columns
-!==========================================================================
-
-!     co2 column at model top (molecule.cm-2)
-
-      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
-                       /(zmmean(lswitch-1)*g)
-
-!     ozone column at model top
-
-      colo3(lswitch-1) = 0.
-
-!     co2 and ozone columns for other levels (molecule.cm-2)
-
-      do l = lswitch-2,1,-1
-         dp = (press(l) - press(l+1))*100.
-         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
-                             *6.022e22*dp/(zmmean(l)*g)
-         col(l) = min(col(l), colairtab(1))
-         colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
-                                 *6.022e22*dp/(zmmean(l)*g)
-      end do
-
-!     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
-
-!     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
-!     profile giving a column 0.1 micron-atmosphere at
-!     a surface pressure of 10 hpa.
-
-      do l = 1,lswitch-1
-         colo3min    = col(l)*7.171e-10
-         ratio_o3(l) = colo3(l)/colo3min
-         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
-         ratio_o3(l) = max(ratio_o3(l), 1.)
-      end do
-
-!==========================================================================
-!     temperature dependence
-!==========================================================================
-
-!     1) search for temperature at 1.9 hPa (tref): vertical interpolation
-
-      tref = temp(1)
-      do l = (lswitch-1)-1,1,-1
-         if (press(l) > 1.9) then
-            cinf = (press(l) - 1.9)/(press(l) - press(l+1))
-            csup = 1. - cinf
-            tref = cinf*temp(l+1) + csup*temp(l)
-            exit
-         end if
-      end do 
-
-!     2) interpolation in temperature
-
-      tref = min(tref,table_temp(1))
-      tref = max(tref,table_temp(ntemp))
-
-      do it = 2, ntemp
-         if (table_temp(it) <= tref) then
-            citemp = (log(tref) - log(table_temp(it)))              &
-                    /(log(table_temp(it-1)) - log(table_temp(it)))
-            indtemp = it - 1
-            exit
-         end if
-      end do
-
-!==========================================================================
-!     loop over vertical levels
-!==========================================================================
-
-      do l = 1,lswitch-1
-
-!     interpolation in air column
-
-         indcol = nz - 1
-         do icol = 1,nz
-            if (colairtab(icol) < col(l)) then
-               indcol = icol - 1
-               exit
-            end if
-         end do
-         cicol = (log(col(l)) - log(colairtab(indcol + 1)))              &
-                /(log(colairtab(indcol)) - log(colairtab(indcol + 1)))
-
-!     interpolation in ozone column
-
-         indozo = nozo - 1
-         do iozo = 1,nozo
-            if (table_ozo(iozo)*10. >= ratio_o3(l)) then
-               indozo = iozo - 1
-               indozo = max(indozo, 1)
-               exit
-            end if
-         end do
-         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)             &
-                /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
-
-!     4-dimensional interpolation weights
-
-!     poids(temp,sza,co2,o3,tau)
-
-         poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau)
-         poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)
-         poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)
-!
-         poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau
-         poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau
-         poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau
-         poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau
-         poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau
-         poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau
-         poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau
-         poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau
-
-!     4-dimensional interpolation in the lookup table
-
-         do ij = 1,nd
-            j(l,ij) =                                                                & 
-            poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij)           &
-          + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)         &
-          + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)         &
-          + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)       &
-          + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)         &
-          + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)       &
-          + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)       &
-          + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)     &
-          + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)         &
-          + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)       &
-          + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)       &
-          + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)     &
-          + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)       &
-          + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)     &
-          + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)     &
-          + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)   &
-!
-          + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)         &
-          + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)       &
-          + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)       &
-          + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)     &
-          + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)       &
-          + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)     &
-          + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)     &
-          + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)   &
-          + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)       &
-          + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)     &
-          + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)     &
-          + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)   &
-          + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)     &
-          + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)   &
-          + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)   &
-          + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
-         end do
-
-!     correction for sun distance
-
-         do ij = 1,nd
-            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
-         end do 
-
-!==========================================================================
-!     end of loop over vertical levels
-!==========================================================================
-
-      end do
-
-      else
-
-!==========================================================================
-!     night
-!==========================================================================
-
-         j(:,:) = 0.
-
-      end if
-
-! photodissociation rates numbering in the lookup table
-
-! jmars.20140930
-
-      j_o2_o         =  1      ! o2 + hv     -> o + o
-      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
-      j_co2_o        =  3      ! co2 + hv    -> co + o
-      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
-      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
-      j_o3_o         =  6      ! o3 + hv     -> o2 + o
-      j_h2o          =  7      ! h2o + hv    -> h + oh
-      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
-      j_ho2          =  9      ! ho2 + hv    -> oh + o
-      j_no           =  10     ! no + hv     -> n + o
-      j_no2          =  11     ! no2 + hv    -> no + o
-      j_hno3         =  12     ! hno3 + hv   -> no2 + oh
-      j_hno4         =  13     ! hno4 + hv   -> no2 + ho2
-
-! jmars.20111014
-
-!     j_o2_o         =  1      ! o2 + hv     -> o + o
-!     j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
-!     j_co2_o        =  3      ! co2 + hv    -> co + o
-!     j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
-!     j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
-!     j_o3_o         =  6      ! o3 + hv     -> o2 + o
-!     j_h2o          =  7      ! h2o + hv    -> h + oh
-!     j_hdo          =  8      ! hdo + hv    -> d + oh
-!     j_h2o2         =  9      ! h2o2 + hv   -> oh + oh
-!     j_ho2          =  10     ! ho2 + hv    -> oh + o
-!     j_no2          =  11     ! no2 + hv    -> no + o
-!     j_ch4_ch3_h    =  12     ! ch4 + hv    -> ch3 + h
-!     j_ch4_1ch2_h2  =  13     ! ch4 + hv    -> 1ch2 + h2
-!     j_ch4_3ch2_h_h =  14     ! ch4 + hv    -> 3ch2 + h + h
-!     j_ch4_ch_h2_h  =  15     ! ch4 + hv    -> ch + h2 + h
-!     j_ch3o2h       =  16     ! ch3o2h + hv -> ch3o + oh
-!     j_ch2o_hco     =  17     ! ch2o + hv   -> h + hco
-!     j_ch2o_co      =  18     ! ch2o + hv   -> h2 + co
-!     j_ch3oh        =  19     ! ch3oh + hv  -> ch3o + h
-!     j_c2h6         =  20     ! c2h6 + hv   -> products
-!     j_hcl          =  21     ! hcl + hv    -> h + cl
-!     j_hocl         =  22     ! hocl + hv   -> oh + cl
-!     j_clo          =  23     ! clo + hv    -> o + cl
-!     j_so2          =  24     ! so2 + hv    -> so + o
-!     j_so           =  25     ! so + hv     -> s + o
-!     j_h2s          =  26     ! h2s + hv    -> hs + s
-!     j_so3          =  27     ! so2 + hv    -> so2 + o
-!     j_hno3         =  28     ! hno3 + hv   -> oh + no2
-!     j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
-
-! fill v_phot array
-
-      v_phot(:,:) = 0.
-
-      do l = 1,lswitch-1
-         v_phot(l, 1) = j(l,j_o2_o)
-         v_phot(l, 2) = j(l,j_o2_o1d)
-         v_phot(l, 3) = j(l,j_co2_o)
-         v_phot(l, 4) = j(l,j_co2_o1d)
-         v_phot(l, 5) = j(l,j_o3_o1d)
-         v_phot(l, 6) = j(l,j_o3_o)
-         v_phot(l, 7) = j(l,j_h2o)
-         v_phot(l, 8) = j(l,j_h2o2)
-         v_phot(l, 9) = j(l,j_ho2)
-         v_phot(l,10) = j(l,j_no)
-         v_phot(l,11) = j(l,j_no2)
-      end do
-
-      return
-      end
-
-!*****************************************************************
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2006)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2007)
@@ -241,6 +241,4 @@
       real zcondicea_co2microp(ngrid,nlayer)
 c     Variables used by the photochemistry
-      logical :: asis             ! true  : adaptative semi-implicit symmetric (asis) chemical solver
-                                  ! false : euler backward chemical solver
       REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
       REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
@@ -1588,16 +1586,4 @@
      $                       surfdust, surfice)
 !           call photochemistry
-!           asis = .false.  euler-backward scheme 
-!           asis = .true.   asis scheme 
-
-            asis = .true.
-
-            if (asis) then
-            call calchim_asis(ngrid,nlayer,nq,
-     &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
-     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
-     $                   zdqcloud,zdqscloud,tauref,co2ice,
-     $                   pu,pdu,pv,pdv,surfdust,surfice)
-            else
             call calchim(ngrid,nlayer,nq,
      &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
