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 "dimensions.h"
!#include "dimphys.h"
#include "chimiedata.h" 
#include "callkeys.h"
!#include "tracer.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 phot(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) 

      return
      end

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*****************************************************************

      implicit none 

!#include "dimensions.h"
!#include "dimphys.h"
#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_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_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    -> cl + o
      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

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.

         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)

      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

      return
      end

c*****************************************************************

      subroutine phot(nlayer,
     $                lswitch, press, temp, sza, tauref, dist_sol, 
     $                rmco2, rmo3, j)

c*****************************************************************

      USE comcstfi_h
      implicit none

!#include "dimensions.h"
!#include "dimphys.h"
#include "chimiedata.h" 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     inputs:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       
      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 :: dist_sol              ! sun distance (AU) 
      real :: rmco2(nlayer)         ! co2 mixing ratio
      real :: rmo3(nlayer)          ! ozone mixing ratio

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     output:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     local:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp, 
     $           iozo, isza, itau, it, l

      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

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     day/night criterion
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      if (sza .le. 95.) then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     temperatures at 1.9 hPa in lookup table
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
      table_temp(1) = 226.2
      table_temp(2) = 206.2
      table_temp(3) = 186.2
      table_temp(4) = 169.8

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     interpolation in solar zenith angle
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      indsza = nsza - 1
      do isza = 1,nsza
         if (szatab(isza) .ge. sza) then
            indsza = min(indsza,isza - 1)
            indsza = max(indsza, 1)
         end if
      end do
      cisza = (sza - szatab(indsza))
     $       /(szatab(indsza + 1) - szatab(indsza))

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     interpolation in dust (tau)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      tau = min(tauref, tautab(ntau))
      tau = max(tau, tautab(1))

      indtau = ntau - 1
      do itau = 1,ntau
         if (tautab(itau) .ge. tau) then
            indtau = min(indtau,itau - 1)
            indtau = max(indtau, 1)
         end if
      end do
      citau = (tau - tautab(indtau))
     $       /(tautab(indtau + 1) - tautab(indtau))

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     co2 and ozone columns
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c     co2 column at model top (molecule.cm-2)

      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.
     $                 /(mugaz*g)

c     ozone column at model top

      colo3(lswitch-1) = 0.

c     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/(mugaz*g)
         col(l) = min(col(l), colairtab(0))
         colo3(l) = colo3(l+1)
     $            + (rmo3(l+1) + rmo3(l))*0.5
     $              *6.022e22*dp/(mugaz*g)
      end do

c     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
c
c     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
c     profile giving a column 0.1 micron-atmosphere at
c     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

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     temperature dependence
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c     1) search for temperature at 1.9 hPa (tref): vertical interpolation

      tref = temp(1)
      do l = (lswitch-1)-1,1,-1
         if (press(l) .gt. 1.9) then
            cinf = (press(l) - 1.9)
     $            /(press(l) - press(l+1))
            csup = 1. - cinf
            tref = cinf*temp(l+1) + csup*temp(l)
            go to 10
         end if
      end do 
 10   continue

c     2) interpolation in temperature

      tref = min(tref,table_temp(1))
      tref = max(tref,table_temp(ntemp))

      do it = 2, ntemp
         if (table_temp(it) .le. tref) then
            citemp = (log(tref) - log(table_temp(it)))
     $              /(log(table_temp(it-1)) - log(table_temp(it)))
            indtemp = it - 1
            goto 20
         end if
      end do
  20  continue

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     loop over vertical levels
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do l = 1,lswitch-1

c     interpolation in air column

         do icol = 0,200
            if (colairtab(icol) .lt. col(l)) then
               cicol = (log(col(l)) - log(colairtab(icol)))
     $                /(log(colairtab(icol-1)) - log(colairtab(icol)))
               indcol = icol - 1
               goto 30
            end if
         end do
 30      continue

cc    interpolation in ozone column

         indozo = nozo - 1
         do iozo = 1,nozo
            if (table_ozo(iozo)*10. .ge. ratio_o3(l)) then
               indozo = min(indozo, iozo - 1)
               indozo = max(indozo, 1)
            end if
         end do
         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)
     $          /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)

cc    4-dimensional interpolation weights

c     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

cc    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

cc    correction for sun distance

         do ij = 1,nd
            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
         end do 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     end of loop over vertical levels
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      end do

      else

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     night
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do ij = 1,nd
         do l = 1,lswitch-1
            j(l,ij) = 0.
         end do
      end do

      end if

      return
      end

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 "dimensions.h"
!#include "dimphys.h"
#include "callkeys.h"
!#include "tracer.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

      return
      end

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 "dimensions.h"
!#include "dimphys.h"
#include "callkeys.h"
!#include "tracer.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

      return
      end

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

!#include "dimensions.h"
!#include "dimphys.h"

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
      return
      end