Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 333)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 334)
@@ -1,5 +1,6 @@
       subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
-     $                   zzlay,zday,pq,pdq,rice,
-     &                   dqchim,dqschim,dqcloud,dqscloud)
+     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
+     $                   dqscloud,tauref,co2ice,
+     $                   pu,pdu,pv,pdv,surfdust,surfice)
 c
       implicit none
@@ -19,4 +20,5 @@
 c    update 03/05/2005 cosmetic changes (Franck Lefevre)
 c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
+c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
 c
 c   Arguments:
@@ -30,9 +32,16 @@
 c    pt(ngridmx,nlayermx)       Temperature (K)
 c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
+c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
+c    pdu(ngridmx,nlayermx)      u component tendency (K)
+c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
+c    pdv(ngridmx,nlayermx)      v component tendency (K)
 c    dist_sol                   distance of the sun (AU)
 c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
 c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
 c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
-c    rice(ngridmx,nlayermx)     Estimated ice crystal radius (m)
+c    tauref(ngridmx)            Optical depth at 7 hPa
+c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
+c    surfdust(ngridmx,nlayermx) dust surface area (micron2/cm3)
+c    surfice(ngridmx,nlayermx)  ice surface area (micron2/cm3)
 c
 c  Output:
@@ -64,12 +73,20 @@
       real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
       real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
+      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
       real    pt(ngridmx,nlayermx)       ! temperature
       real    pdt(ngridmx,nlayermx)      ! temperature tendency
+      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
+      real    pdu(ngridmx,nlayermx)      ! u component tendency
+      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
+      real    pdv(ngridmx,nlayermx)      ! v component tendency
       real    dist_sol                   ! distance of the sun (AU)
-      real    mu0(ngridmx)       ! cos of solar zenith angle (=1 when sun at zenith)
+      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
       real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
       real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
       real    zday                       ! date (time since Ls=0, in martian days)
-      real    rice(ngridmx,nlayermx)     ! Estimated ice crystal radius (m)
+      real    tauref(ngridmx)            ! optical depth at 7 hPa
+      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
+      real    surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3)
+      real    surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3)
 
 c   outputs:
@@ -87,8 +104,15 @@
       integer  ig,l,i,iq
       integer  foundswitch, lswitch
+      integer  ig_vl1
+      real latvl1, lonvl1
       real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
                                      ! new mole fraction after
-      real colden(ngridmx,nqmx)      ! Column densities (cm-2)
       real zt(ngridmx,nlayermx)      ! temperature
+      real zu(ngridmx,nlayermx)      ! u component of the wind
+      real zv(ngridmx,nlayermx)      ! v component of the wind
+      real taucol                    ! optical depth at 7 hPa
+
+      logical depos                  ! switch for dry deposition
+      parameter (depos=.false.)
 c
 c  for each column of atmosphere:
@@ -100,11 +124,11 @@
       real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
       real szacol                 !  Solar zenith angle
+      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
+      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
       real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
 c
 c for output:
 c
-      real zdens3d(ngridmx,nlayermx) !  Density  (cm-3)
       real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
-      real surfice(ngridmx,nlayermx) !  Surface of ice particules (um2/cm3)
 
       logical output              ! to issue calls to writediagfi and stats
@@ -112,8 +136,9 @@
 
       logical,save :: firstcall=.true.
-      integer,save :: nbq  ! number of tracers used in the chemistry
+      integer,save :: nbq       ! number of tracers used in the chemistry
       integer,save :: niq(nqmx) ! array storing the indexes of the tracers
 
 ! index of tracers:
+
       integer,save :: i_co2=0
       integer,save :: i_co=0
@@ -127,13 +152,9 @@
       integer,save :: i_ho2=0
       integer,save :: i_h2o2=0
+      integer,save :: i_ch4=0
       integer,save :: i_n2=0
       integer,save :: i_ar=0
       integer,save :: i_ice=0 ! water ice
       integer,save :: i_h2o=0 ! water vapour
-
-c
-c  scheme A: 1 ; scheme B: 2
-c
-      integer,parameter :: scheme=2
 c
 c=======================================================================
@@ -144,13 +165,6 @@
 c
          if (photochem) then
-            print*,'calchim: INIT CHEMISTRY'
-            if (scheme  .eq.  1) then
-               print*,'calchim: Scheme A : A METTRE A JOUR !!'
-               stop
-c              call init_chimie_A
-            else
-               print*,'calchim: Scheme B'
-               call init_chimie_B
-            end if
+            print*,'calchim: Read photolysis lookup table'
+            call read_phototable
          end if
 
@@ -245,4 +259,12 @@
            niq(nbq)=i_h2o2
          endif
+         i_ch4=igcm_ch4
+         if (i_ch4.eq.0) then
+           write(*,*) "calchim: Error; no CH4 tracer !!!"
+           stop
+         else
+           nbq=nbq+1
+           niq(nbq)=i_ch4
+         endif
          i_n2=igcm_n2
          if (i_n2.eq.0) then
@@ -278,22 +300,21 @@
          endif
 
-         write(*,*) 'calchim: found nbq=',nbq,' tracers'
-         write(*,*) '         i_co2=',i_co2
-         write(*,*) '         i_co=',i_co
-         write(*,*) '         i_o=',i_o
-         write(*,*) '         i_o1d=',i_o1d
-         write(*,*) '         i_o2=',i_o2
-         write(*,*) '         i_o3=',i_o3
-         write(*,*) '         i_h=',i_h
-         write(*,*) '         i_h2=',i_h2
-         write(*,*) '         i_oh=',i_oh
-         write(*,*) '         i_ho2=',i_ho2
-         write(*,*) '         i_h2o2=',i_h2o2
-         write(*,*) '         i_n2=',i_n2
-         write(*,*) '         i_ar=',i_ar
-         write(*,*) '         i_ice=',i_ice
-         write(*,*) '         i_h2o=',i_h2o
-!         write(*,*) '         niq(:)=',niq
-!         write(*,*) '     nqchem_min,nqmx=',nqchem_min,nqmx
+         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
+         write(*,*) '               i_co2  = ',i_co2
+         write(*,*) '               i_co   = ',i_co
+         write(*,*) '               i_o    = ',i_o
+         write(*,*) '               i_o1d  = ',i_o1d
+         write(*,*) '               i_o2   = ',i_o2
+         write(*,*) '               i_o3   = ',i_o3
+         write(*,*) '               i_h    = ',i_h
+         write(*,*) '               i_h2   = ',i_h2
+         write(*,*) '               i_oh   = ',i_oh
+         write(*,*) '               i_ho2  = ',i_ho2
+         write(*,*) '               i_h2o2 = ',i_h2o2
+         write(*,*) '               i_ch4  = ',i_ch4
+         write(*,*) '               i_n2   = ',i_n2
+         write(*,*) '               i_ar   = ',i_ar
+         write(*,*) '               i_ice  = ',i_ice
+         write(*,*) '               i_h2o  = ',i_h2o
          
          firstcall = .false.
@@ -302,8 +323,12 @@
 ! Initialize output tendencies to zero (to handle case of tracers which
 ! are not used in the chemistry (e.g. dust))
-      dqchim(:,:,:)=0
-      dqschim(:,:)=0
-
-
+
+      dqchim(:,:,:) = 0
+      dqschim(:,:)  = 0
+c
+!     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.)
 c
 c=======================================================================
@@ -317,10 +342,13 @@
          foundswitch = 0
          do l = 1,nlayermx
-            zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
-            do i=1,nbq
-              iq=niq(i) ! get tracer index
+            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)
-            enddo
+              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
+c
             zpress(l) = pplay(ig,l)/100.
             ztemp(l)  = zt(ig,l)
@@ -328,5 +356,10 @@
             zlocal(l) = zzlay(ig,l)/1000.
 c
-c     search for switch index between regions
+c           surfdust1d and surfice1d in cm2/cm3
+c
+            surfdust1d(l) = surfdust(ig,l)*1.e-8
+            surfice1d(l)  = surfice(ig,l)*1.e-8
+c
+c           search for switch index between regions
 c
             if (photochem .and. thermochem) then
@@ -340,24 +373,11 @@
             end if
             if (.not.  thermochem) then
-               lswitch = min(33,nlayermx+1)
+               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
             end if
 c
-c     ice surface area  in microns^2/cm^3
-c
-c     = 4 pi r^2 * [ zq * mugaz/NA / (rhoice*4/3 pi r^3) ] *zdens
-c     = 3/r * [ zq * mugaz/NA / rhoice ] *zdens
-c     with r in microns, rhoice = 0.92e-12 g microns^-3 and zdens in cm^-3
-c
-            if (water) then
-               zycol(l,i_ice) = (3.e-6/rice(ig,l))*zq(ig,l,i_ice)
-     $                           *(mugaz/6.022e23)*zdens(l)/0.92e-12
-c              write(*,*) "rice=",rice(ig,l)," m / zdens=",zdens(l),
-c    $        " cm-3 / icesurf=",zycol(l,nqmx-1)," microns^2/cm^3"
-               surfice(ig,l) = zycol(l,i_ice)
-            end if
-c
          end do ! of do l=1,nlayermx
 c
          szacol = acos(mu0(ig))*180./pi
+         taucol = tauref(ig)
 c
 c=======================================================================
@@ -365,19 +385,25 @@
 c=======================================================================
 c
-        if (photochem) then
-           if (scheme .eq. 1) then
-              print*,'Scheme A : A METTRE A JOUR !!'
-c             call photochemist_A(zycol,szacol,ptimestep,
-c    $                            zpress,ztemp,zdens,dist_sol)
-           else
-              call photochemist_B(lswitch,zycol,szacol,ptimestep,
-     $                            zpress,ztemp,zdens,dist_sol,jo3)
-           end if
-        end if
-
-        if (thermochem) then
-           call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
-     $                      zlocal,szacol,ptimestep,zday)
-        end if
+c        chemistry in lower atmosphere
+c
+         if (photochem) then
+            call photochemistry(lswitch,zycol,szacol,ptimestep,
+     $                          zpress,ztemp,zdens,dist_sol,
+     $                          surfdust1d,surfice1d,jo3,taucol)
+         end if
+c
+c        chemistry in upper atmosphere
+c
+         if (thermochem) then
+            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
+     $                       zlocal,szacol,ptimestep,zday)
+         end if
+c
+c        dry deposition
+c
+         if (depos) then
+            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev, 
+     $                      zu, zv, zt, zycol, ptimestep, co2ice)
+         end if
 c
 c=======================================================================
@@ -389,5 +415,4 @@
          if (water) then
             do l = 1,nlayermx
-!               dqchim(ig,l,nqmx-1) = 0.
                dqchim(ig,l,i_ice) = 0.
             end do
@@ -425,5 +450,4 @@
      $                                 - zq(ig,l,iq))/ptimestep
                    else if (iq.eq.i_o) then
-!                      i_o = iq
                       dqchim(ig,l,iq) = 0.
                    end if
@@ -440,15 +464,4 @@
           end do ! of do l = 1,nlayermx
 c
-c     dust: This is now taken care of as a first step at beginning of routine
-c
-!          if (nqchem_min .gt. 1) then
-!             do iq = 1,nqchem_min-1
-!                do l = 1,nlayermx
-!                   dqchim(ig,l,iq) = 0.
-!                end do
-!                dqschim(ig,iq) = 0.
-!             end do
-!          end if
-c
 c     condensation of h2o2
 c
@@ -456,44 +469,7 @@
      $                 ztemp,zycol,dqcloud,dqscloud)
 c
-c     for outputs
-c
-          do i=1,nbq
-            iq=niq(i) ! get tracer index
-             colden(ig,iq) = 0.
-             do l = 1,nlayermx
-c
-c     column density converted in cm-2
-c     pplev en pa, mugaz en g.mol-1 et g en m.s-2
-c     not for ice
-c
-                if (iq.ne.i_h2o2) then
-                  colden(ig,iq) = colden(ig,iq) + zycol(l,iq)
-     $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
-     $                         /(mmean(ig,l)*g)
-                else   ! for H2O2, remove condensation from zycol
-                  colden(ig,iq) = colden(ig,iq) + (zycol(l,iq) +
-     $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
-     $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
-     $                         /(mmean(ig,l)*g)
-                end if
-c
-c     local densities, for outputs (put in zq)
-c     not for ice
-c
-                zq(ig,l,iq) = zycol(l,iq)*zdens(l)
-c                        for H2O2, remove condensation from zycol
-                if (iq.eq.i_h2o2) then
-                   zq(ig,l,iq) = zdens(l)*(zycol(l,iq) +
-     $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
-                end if
-             end do
-          end do
-c
 c     density and j(o3->o1d), for outputs
 c
-          zdens3d(ig,1) = zdens(1)
-          jo3_3d(ig,1) = jo3(1)
-          do l = 2,nlayermx
-             zdens3d(ig,l) = zdens(l)
+          do l = 1,nlayermx
              jo3_3d(ig,l) = jo3(l)
           end do
@@ -513,50 +489,10 @@
 
       if (output) then
-
          if (ngridmx .gt. 1) then
-c           call writediagfi(ngridmx,'dens','atm dens.','cm-3',3,zdens3d(1,1))
-c           call writediagfi(ngridmx,'jo3','j o3->o1d','s-1',3,jo3_3d(1,1))
-c           call writediagfi(ngridmx,'sice','ice surf.','um2/cm3',3,surfice(1,1))
-            do i=1,nbq
-              iq=niq(i) ! get tracer index
-               if (iq.ne.i_ice) then
-                 write(str20(1:20),'(a20)') noms(iq)
-                 call writediagfi(ngridmx,'n_'//trim(str20),'density',
-     $                             'cm-3',3,zq(1,1,iq))
-c                 call writediagfi(ngridmx,'dqch_'//str5,'density','cm-3',3,dqchim(1,1,iq))
-c                 if (noms(iq) .eq. "h2o2" .or. noms(iq) .eq. "h2o") then
-c                    call writediagfi(ngridmx,'cl_'//str5,'density','cm-3',3,dqcloud(1,1,iq))
-c                 end if
-                 call writediagfi(ngridmx,'c_'//trim(str20),
-     $                            'col. dens.','cm-2',2,colden(1,iq))
-               end if
-            end do
-c
-            if (callstats) then
-c
-c              convert to mole.cm-2 for the column densities
-c
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  do ig = 1,ngridmx
-                     colden(ig,iq) = colden(ig,iq)/6.022e23
-                  end do   
-               end do 
-c
-c              call wstats(ngridmx,"jo3","jo3->o1d","s-1",3,jo3_3d)
-c
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  if (iq.ne.i_ice) then
-                     write(str20(1:20),'(a20)') noms(iq)
-                     call wstats(ngridmx,"n_"//trim(str20),"density",
-     &                           "cm-3",3,zq(1,1,iq))
-                     call wstats(ngridmx,"c_"//trim(str20),"col. dens.",
-     &                           "mol cm-2",2,colden(1,iq))
-                  end if
-               end do ! of i=1,nbq
-            end if ! of if (callstats)
+            call writediagfi(ngridmx,'jo3','j o3->o1d',
+     $                       's-1',3,jo3_3d(1,1))
+            call wstats(ngridmx,'jo3','j o3->o1d',
+     $                       's-1',3,jo3_3d(1,1))
          end if ! of if (ngridmx.gt.1)
-c
       endif ! of if (output)
 c
Index: trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h	(revision 333)
+++ trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h	(revision 334)
@@ -2,59 +2,40 @@
 c ---  Data for chemical routines  ---
 c ------------------------------------
-c    Author: Sebastien Lebonnois
-c            Franck Lefevre
-c    
-c  Update after water clouds included
-c ------------------------------------
 c
-      integer    nd, nozo, nr, nsza, ntemp
+      integer nd, nozo, nr, nsza, ntemp, ntau
 c
-      parameter (nd    = 11)
+      parameter (nd    = 29)
       parameter (nozo  = 7)
       parameter (nr    = nd + 28)
       parameter (nsza  = 27)
       parameter (ntemp = 4)
+      parameter (ntau  = 8)
 c
-      real       kb
-      parameter (kb = 1.38e-23)
+      real kb
+      parameter (kb = 1.3806e-23)
 c
-c     factor: molecular gas constant/(boltzmann*gravity)
+      common/chimiedata/jphot,colairtab,table_ozo
 c
-      real       factor
-      parameter (factor = 192./(kb*1.e4*3.72))
+      real jphot(ntemp,nsza,0:200,nozo,ntau,nd)
+      real colairtab(0:200)
+      real szatab(nsza)
+      real table_ozo(nozo)
+      real tautab(ntau)
 c
-      common/chimiedata/reactifs,nprod,nperte,prod,perte,jphot,
-     $                  colairtab,table_ozo
+      data szatab/0.,  5., 10., 15., 20., 25.,
+     $            30., 35., 40., 45., 50., 55.,
+     $            60., 65., 70., 75., 80., 82.,
+     $            84., 86., 88., 90., 91., 92.,
+     $            93., 94., 95./
 c
-c     used in obsolete scheme A :
-      integer    reactifs(nr,5),nprod(nqmx),nperte(nqmx)
-      integer    prod(nqmx,100),perte(nqmx,100,2)
-c
-c     VERSION: 28/05/2004  (update: Ar)
-c
-c     photodissociation rates matrix: version dependent
-c     for this version, file jmars.20030707
-c
-      real       jphot(ntemp,nsza,0:200,nozo,nd)
-      real       colairtab(0:200)
-      real       szatab(nsza)
-      real       table_ozo(nozo)
-c
-      data       szatab/0.,  5., 10., 15., 20., 25.,
-     $                 30., 35., 40., 45., 50., 55.,
-     $                 60., 65., 70., 75., 80., 82.,
-     $                 84., 86., 88., 90., 91., 92.,
-     $                 93., 94., 95./
-c
-c     compounds taken into account: version dependent
-c     DOES NOT INCLUDE h2o NOR water ice
+      data tautab/0., 0.2, 0.4, 0.6, 0.8, 1., 2., 4./
 c
       integer       ncomp
-      parameter    (ncomp=13)
+      parameter    (ncomp = 15)
       character*10  nomchem(ncomp)
       data          nomchem/"co2","co","o","o1d","o2","o3","h","h2",
-     $                   "oh","ho2","h2o2", "n2", "ar"/
+     $                      "oh","ho2","h2o2", "ch4", "n2", "ar", "h2o"/
 c
-      real  mmolchem(ncomp)
-      data  mmolchem/44.,28.,16.,16.,32.,48.,1.,2.,17.,33.,34.,28.,40./
-
+      real mmolchem(ncomp)
+      data mmolchem/44.,28.,16.,16.,32.,48., 1.,
+     $               2.,17.,33.,34.,16.,28.,40.,18./
Index: trunk/LMDZ.MARS/libf/aeronomars/init_chimie_B.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/init_chimie_B.F	(revision 333)
+++ 	(revision )
@@ -1,83 +1,0 @@
-      subroutine init_chimie_B
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c   Scheme B
-c
-c   subject:
-c   --------
-c
-c   Initialization of the photochemical scheme
-c
-c   VERSION: 30/04/2003
-c
-c   This is dependent of the chemistry version:
-c   for scheme B its only purpose is reading the
-c   photodissociation rates lookup table
-c
-c   Author:   Franck Lefevre for scheme B
-c
-c   Arguments:
-c   ----------
-c
-c   The output variable is jphot and is put in common chimiedata.
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "chimiedata.h"
-#include "datafile.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer fic, ij, iozo, isza, itemp, iz, ierr
-c
-      real xsza
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      fic = 81
-
-      open(fic, form = 'formatted', status = 'old',
-     &file =trim(datafile)//'/jmars.20030707',iostat=ierr)
-
-      IF (ierr.NE.0) THEN
-        write(*,*)'Error : cannot open file jmars.20030707 '
-        write(*,*)'(in aeronomars/init_chimie_B.F)'
-        write(*,*)'It should be in :',trim(datafile),'/'
-        write(*,*)'1) You can change this directory address in '
-        write(*,*)'   file phymars/datafile.h'
-        write(*,*)'2) If necessary jmars.20030707 (and other datafiles)'
-        write(*,*)'   can be obtained online on:'
-        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
-        STOP
-      ENDIF
-
-c
-      print*, 'init_chimie_B: lecture table des J...'
-c
-      do itemp = 1,ntemp
-         do iozo = 1,nozo
-            do isza = 1,nsza
-               do iz = 200,0,-1
-                  read(fic,*) colairtab(iz), xsza, table_ozo(iozo)
-                  read(fic,*) (jphot(itemp,isza,iz,iozo,ij), ij=1,7)
-                  read(fic,*) (jphot(itemp,isza,iz,iozo,ij), ij=8,11)
-                  do ij = 1,nd
-                     if (jphot(itemp,isza,iz,iozo,ij) .eq. 1.e-30) then
-                        jphot(itemp,isza,iz,iozo,ij) = 0.
-                     end if
-                  end do
-               end do
-            end do
-         end do
-      end do
-c
-      print*, 'init_chimie_B: lecture table des J...ok'
-      close(fic)
-c
-      return
-      end
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemist_B.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemist_B.F	(revision 333)
+++ 	(revision )
@@ -1,1710 +1,0 @@
-c*****************************************************************
-c
-c     Photochemical routines for Scheme B
-c
-c     Author: Franck Lefevre
-c     ------
-c
-c     update 12/06/2003: compatibility with dust and iceparty
-c        (S. Lebonnois)
-c     update sept. 2008: identify tracer indexes using tracer names
-c        (Ehouarn Millour)
-c
-c*****************************************************************
-c
-      subroutine photochemist_B(lswitch,zycol, sza, ptimestep, press,
-c     $                          temp, dens, dist_sol)
-     $                          temp, dens, dist_sol, jo3)
-c
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "chimiedata.h" 
-#include "callkeys.h"
-#include "tracer.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     input/output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real zycol(nlayermx,nqmx)  ! chemical species volume mixing ratio
-c if water, zycol(l,i_h2o_ice) is ice surface area (in microns^2/cm^3)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real sza                   ! solar zenith angle (deg)
-      real ptimestep             ! physics timestep (s)
-      real press(nlayermx)       ! pressure (hpa)
-      real temp(nlayermx)        ! temperature (k)
-      real dens(nlayermx)        ! density (cm-3)
-      real dist_sol              ! sun distance (au) 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real jo3(nlayermx)  ! photodissociation rate O3 -> O1D
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer phychemrat         ! ratio physics timestep/chemistry timestep
-      integer istep
-      integer i_o3,lev,iq
-      integer nesp
-      integer lswitch
-c
-      parameter (nesp = 15)      ! number of species in the chemistry
-c
-      real ctimestep             ! chemistry timestep (s) 
-      real rm(nlayermx,nesp)     ! species volume mixing ratio 
-      real j(nlayermx,nd)        ! interpolated photolysis rates (s-1)
-      real icesurf(nlayermx)     ! ice surface area (cm^2/cm^3)
-      real rmo3(nlayermx)        ! ozone mixing ratio
-c
-c     reaction rates
-c
-      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
-      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
-     $     b004(nlayermx), b005(nlayermx), b006(nlayermx)
-      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
-     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
-     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
-     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
-     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
-     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
-      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
-      real e001(nlayermx), e002(nlayermx)
-      real h001(nlayermx), h002(nlayermx)
-
-      logical,save :: firstcall=.true.
-      integer,save :: i_h2o_ice=0 ! water ice tracer index
-
-! Initializations (at first call only)
-      if (firstcall) then
-        ! find index of water ice tracer
-        do iq=1,nqmx
-          if (noms(iq).eq."h2o_ice") then
-            i_h2o_ice=iq
-          endif
-        enddo
-        ! verify that it is the same index as the global one:
-        if (i_h2o_ice.ne.igcm_h2o_ice) then
-          write(*,*) "photochemist_B: Error ; found water ice tracer",
-     &               " index ",i_h2o_ice," does not match GCM index ",
-     &               " igcm_h2o_ice=",igcm_h2o_ice
-          stop
-        endif
-        ! also check that there indeed is an "h2o_ice" tracer
-        if (i_h2o_ice.eq.0) then
-          write(*,*) "photochemist_B: Error ; did not find",
-     &               " water ice tracer"
-          stop
-        endif
-        firstcall=.false.
-      endif
-
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     ctimestep : chemistry timestep (s)                         c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      ctimestep = 10.*60.
-c
-      phychemrat = int(ptimestep/ctimestep)
-c     print*,"PHYCHEMRAT = ",phychemrat
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialisation of chemical species and families            c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call gcmtochim(zycol, lswitch, nesp, rm) 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     ice surface area for heterogenous chemistry
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      if (water) then
-       do lev=1,lswitch-1
-!        icesurf(lev)= zycol(lev,nqmx-1)*1.e-8 ! microns^2 -> cm^2
-        icesurf(lev)= zycol(lev,i_h2o_ice)*1.e-8 ! microns^2 -> cm^2
-       enddo
-      else
-       do lev=1,lswitch-1
-        icesurf(lev)=0.0e0
-       enddo
-      endif
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     compute reaction rates                                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call chemrates(lswitch, dens, press, temp, icesurf,
-     $               a001, a002, a003,
-     $               b001, b002, b003, b004, b005, b006,
-     $               c001, c002, c003, c004, c005, c006,
-     $               c007, c008, c009, c010, c011, c012,
-     $               c013, c014, c015, c016, c017, c018,
-     $               d001, d002, d003,
-     $               e001, e002,
-     $               h001, h002)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     interpolation of photolysis rates in the lookup table      c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      i_o3 = 5 
-c
-      do lev=1,lswitch-1
-        rmo3(lev)=rm(lev,i_o3)
-      enddo
-
-      call phot(lswitch, press, temp, sza, dist_sol, rmo3, j)
-c
-      do lev=1,lswitch-1
-        jo3(lev) = j(lev,5)
-      enddo
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over time within the physical timestep                c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do istep = 1,phychemrat
-         call chimie_B(lswitch,nesp, rm, j, dens, ctimestep, 
-     $                 press, temp, sza, dist_sol,
-     $                 a001, a002, a003,
-     $                 b001, b002, b003, b004, b005, b006,
-     $                 c001, c002, c003, c004, c005, c006,
-     $                 c007, c008, c009, c010, c011, c012,
-     $                 c013, c014, c015, c016, c017, c018,
-     $                 d001, d002, d003,
-     $                 e001, e002,
-     $                 h001, h002)
-c
-c        print*,'appel chimie ok iteration ', istep
-c
-      end do
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     save chemical species for the gcm                          c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call chimtogcm(zycol, lswitch, nesp, rm) 
-c
-      return
-      end
-c
-c*****************************************************************
-c
-      subroutine chimie_B(lswitch, nesp, rm, j, dens, dt, 
-     $                    press, t, sza, dist_sol,
-     $                    a001, a002, a003,
-     $                    b001, b002, b003, b004, b005, b006,
-     $                    c001, c002, c003, c004, c005, c006,
-     $                    c007, c008, c009, c010, c011, c012,
-     $                    c013, c014, c015, c016, c017, c018,
-     $                    d001, d002, d003,
-     $                    e001, e002,
-     $                    h001, h002)
-c
-c*****************************************************************
-c
-      implicit none 
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "chimiedata.h" 
-#include "callkeys.h" 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     input/output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      integer lswitch         ! interface level between chemistries
-      integer nesp            ! number of species
-c
-      real rm(nlayermx,nesp)  ! volume mixing ratios 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real dens(nlayermx)     ! density (cm-3)
-      real dt                 ! chemistry timestep (s) 
-      real j(nlayermx,nd)     ! interpolated photolysis rates (s-1) 
-      real press(nlayermx)    ! pressure (hpa)
-      real t(nlayermx)        ! temperature (k)
-      real sza                ! solar zenith angle (deg)
-      real dist_sol           ! sun distance (au) 
-c
-c     reaction rates
-c
-      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
-      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
-     $     b004(nlayermx), b005(nlayermx), b006(nlayermx)
-      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
-     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
-     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
-     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
-     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
-     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
-      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
-      real e001(nlayermx), e002(nlayermx)
-      real h001(nlayermx), h002(nlayermx)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer HETERO
-      parameter (HETERO=0)
-
-      integer iesp, iter, l, niter
-      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_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
-c
-      parameter (niter = 5)          ! iterations in the chemical scheme
-c
-      real*8 cc0(nlayermx,nesp)        ! initial number density (cm-3)
-      real*8 cc(nlayermx,nesp)         ! final number density (cm-3)
-      real*8 co2(nlayermx)             ! co2 number density (cm-3)
-      real*8 nox(nlayermx)             ! nox number density (cm-3)
-      real*8 no(nlayermx)              ! no  number density (cm-3)
-      real*8 no2(nlayermx)             ! no2 number density (cm-3)
-      real*8 production(nlayermx,nesp) ! production rate
-      real*8 loss(nlayermx,nesp)       ! loss rate
-c
-      real c007h(nlayermx)           ! c007 rate modified
-c                                    ! by heterogenous reactions h001 and h002
-c
-      real ro_o3, rh_ho2, roh_ho2, rno2_no
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     tracer numbering in the chemistry
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      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_h2o    = 12
-      i_n2     = 13
-      i_hox    = 14
-      i_ox     = 15
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     numbering of photolysis rates
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      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
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     volume mixing ratio -> density conversion
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      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
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     co2 and nox number densities (cm-3)    
-c      
-c     co2 mixing ratio: 0.953
-c     nox mixing ratio: 6.e-10 (yung and demore, 1999)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-         co2(l) = cc(l,i_co2)
-         nox(l) = 6.e-10*dens(l)
-      end do 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over iterations in the chemical scheme
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do iter = 1,niter
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     modification of c007
-c        by heterogenous reactions h001 and h002
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      if (HETERO.eq.1) then
-       do l = 1,lswitch-1
-         c007h(l) = c007(l)
-     $            + h001(l)/max(cc(l,i_oh), 1.d-30*dens(l))
-     $            + h002(l)/max(cc(l,i_ho2),1.d-30*dens(l))
-c         write(*,*) "C007 = ",c007(l)," / C007H = ",c007h(l)
-       end do
-      else
-       do l = 1,lswitch-1
-         c007h(l) = c007(l)
-       end do
-      endif
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        nox species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no2/no
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-c
-         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.d-30*dens(l)))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         no(l) = nox(l)/(1. + rno2_no)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        no2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         no2(l) = no(l)*rno2_no
-c
-      end do
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        hox species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        photochemical equilibrium for oh and ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h/ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-c
-         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)
-     $           + c007h(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)
-     $           + j(l,j_ho2))
-     $          /(c011(l)*cc(l,i_o2))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oh/ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         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)
-     $            + j(l,j_h2o)*cc(l,i_h2o)
-     $             /max(cc(l,i_ho2),dens(l)*1.d-30))
-     $            /(c002(l)*cc(l,i_o)
-     $            + c007h(l)*cc(l,i_ho2)
-     $            + e001(l)*cc(l,i_co))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         cc(l,i_h) = cc(l,i_hox)
-     $                 /(1. + (1. + roh_ho2)/rh_ho2)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         cc(l,i_ho2) = cc(l,i_h)/rh_ho2
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oh
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
-c
-      end do
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     modification of c007
-c        by heterogenous reactions h001 and h002
-c           after update of cc(l,i_oh) and cc(l,i_ho2)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      if (HETERO.eq.1) then
-       do l = 1,lswitch-1
-         c007h(l) = c007(l)
-     $            + h001(l)/max(cc(l,i_oh), 1.d-30*dens(l))
-     $            + h002(l)/max(cc(l,i_ho2),1.d-30*dens(l))
-c         write(*,*) "C007 = ",c007(l)," / C007H = ",c007h(l)
-       end do
-      endif
-c
-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
-c        
-      if (sza .le. 95.) then
-c
-         do l = 1,lswitch-1
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o(1d)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            cc(l,i_o1d) = (j(l,j_co2_o1d)*co2(l)
-     $                   + j(l,j_o2_o1d)*cc(l,i_o2)
-     $                   + j(l,j_o3_o1d)*cc(l,i_o3))
-     $                   /( b001(l)*co2(l)
-     $                   + 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))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o/o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            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))
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            cc(l,i_o) = cc(l,i_o3)*ro_o3
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        ox = o + o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            production(l,i_ox) =
-     $                   + j(l,j_co2_o)*co2(l)
-     $                   + j(l,j_co2_o1d)*co2(l)
-     $                   + 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)
-c
-            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)
-c
-            loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
-c
-         end do
-c
-      else
-c
-         do l = 1,lswitch-1
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o(1d)
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            cc(l,i_o1d) = 0.
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o3
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            production(l,i_o3) = a001(l)*cc(l,i_o2)*cc(l,i_o)
-c
-            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)
-c
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-            production(l,i_o) = c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                        + c013(l)*cc(l,i_oh)*cc(l,i_oh)
-c
-            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)
-c
-         end do
-      end if
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     other species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        co
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         production(l,i_co) = j(l,j_co2_o)*co2(l)
-     $                      + j(l,j_co2_o1d)*co2(l)
-c
-         loss(l,i_co) = e001(l)*cc(l,i_oh)
-     $                + e002(l)*cc(l,i_o)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        o2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         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)
-     $                + c007h(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)
-c
-         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)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         production(l,i_h2) = c005(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                      + c018(l)*cc(l,i_h)*cc(l,i_h)
-c
-         loss(l,i_h2) = b003(l)*cc(l,i_o1d)
-     $                + c010(l)*cc(l,i_oh)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2o
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         production(l,i_h2o) = 
-     $                   c006(l)*cc(l,i_h)*cc(l,i_ho2)
-     $                 + c007h(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)
-c
-         loss(l,i_h2o) = j(l,j_h2o)
-     $                 + b002(l)*cc(l,i_o1d)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        h2o2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         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
-         loss(l,i_h2o2) = j(l,j_h2o2)
-     $                  + c009(l)*cc(l,i_oh)
-     $                  + c012(l)*cc(l,i_o)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        hox = h + oh + ho2
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-         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)
-c
-         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.*c007h(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)
-c
-         loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
-c
-      end do
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     update number densities
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c     long-lived species
-c
-      do l = 1,lswitch-1
-         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)
-      end do
-c
-c     ox species
-c
-      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 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     end of loop over iterations
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      end do
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     density -> volume mixing ratio conversion
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do iesp = 1,nesp
-         do l = 1,lswitch-1
-            rm(l,iesp) = max(cc(l,iesp)/dens(l), 1.d-30)
-         end do
-      end do
-c
-      return
-      end
-c
-c*****************************************************************
-c
-      subroutine phot(lswitch, press, temp, sza, dist_sol, rmo3, j)
-c
-c*****************************************************************
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "chimiedata.h" 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c       
-      integer lswitch         ! interface level between chemistries
-      real press(nlayermx)             ! pressure (hPa)
-      real temp(nlayermx)              ! temperature (K)
-      real sza                         ! solar zenith angle (deg)
-      real dist_sol                    ! sun distance (AU) 
-      real rmo3(nlayermx)              ! ozone mixing ratio
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real j(nlayermx,nd)              ! interpolated photolysis rates (s-1) 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer icol, ij, indsza, indcol, indozo, indtemp, 
-     $        iozo, isza, it, iztop, l
-c
-      real col(nlayermx)               ! overhead air column   (molecule cm-2)
-      real colo3(nlayermx)             ! overhead ozone column (molecule cm-2)
-      real poids(2,2,2,2)              ! 4D 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
-      real colo3min, dp, scaleh
-      real ratio_o3(nlayermx)
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     day/night criterion
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      if (sza .le. 95.) then 
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     temperatures at 1.9 hPa in jmars
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      table_temp(1) = 226.2
-      table_temp(2) = 206.2
-      table_temp(3) = 186.2
-      table_temp(4) = 169.8
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     interpolation in solar zenith angle
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      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))
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     air and ozone columns
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c     air column at model top: assign standard value
-c
-      scaleh = 10.
-      iztop = nint(scaleh*log(press(1)/press(lswitch-1)))
-      iztop = min(iztop,200)
-      col(lswitch-1) = colairtab(iztop)
-c
-c     ozone column at model top
-c
-      colo3(lswitch-1) = 0.
-c
-c     air and ozone columns for other levels
-c
-      do l = lswitch-2,1,-1
-         dp = (press(l) - press(l+1))*100.
-         col(l) = col(l+1) + factor*dp
-         col(l) = min(col(l), colairtab(0))
-         colo3(l) = colo3(l+1)
-     $            + factor*(rmo3(l+1) + rmo3(l))*0.5*dp
-      end do
-c
-c     ratio ozone column/minimal theoretical column (ro3 = 7.171e-10)
-c
-      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), table_ozo(1)*10.)
-      end do
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     temperature dependence
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c     1) search for temperature at 1.9 hPa (tref): vertical interpolation
-c
-c      print*,'****** recherche de T a 1.9 hPa ******'
-      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)
-c
-c            print*, press(l+1), temp(l+1)
-c            print*, press(l  ), temp(l  )
-c            print*, 'tref = ', tref
-            go to 10
-         end if
-      end do 
- 10   continue
-c
-c     2) interpolation in temperature
-c
-      tref = min(tref,table_temp(1))
-      tref = max(tref,table_temp(ntemp))
-c
-      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
-c            print*,'interpolation entre ', table_temp(it), ' et ',
-c     $              table_temp(it-1) 
-c            print*,'indtemp = ', indtemp, ' citemp = ', citemp
-            goto 20
-         end if
-      end do
-  20  continue
-c      print*,'****** fin interpolation en T ********'
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over vertical levels
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-c
-c     interpolation in air column
-c
-         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
-c
-cc    interpolation in ozone column
-c
-         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.)
-c
-c         write(6,'(i2,4e13.4,5f10.4)') 
-c    $         l, press(l), col(l), 
-c    $         colairtab(indcol), colairtab(indcol+1), 
-c    $         colo3(l)/2.687e15, 
-c    $         ratio_o3(l), table_ozo(indozo), table_ozo(indozo+1), 
-c    $         ciozo
-c 
-cc    4-dimensional interpolation weights
-c
-         poids(1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)
-         poids(1,1,1,2) = citemp*(1.-cisza)*cicol*ciozo
-         poids(1,1,2,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)
-         poids(1,1,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo
-         poids(1,2,1,1) = citemp*cisza*cicol*(1.-ciozo)
-         poids(1,2,1,2) = citemp*cisza*cicol*ciozo
-         poids(1,2,2,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)
-         poids(1,2,2,2) = citemp*cisza*(1.-cicol)*ciozo
-         poids(2,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)
-         poids(2,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo
-         poids(2,1,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)
-         poids(2,1,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo
-         poids(2,2,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)
-         poids(2,2,1,2) = (1.-citemp)*cisza*cicol*ciozo
-         poids(2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)
-         poids(2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo
-c
-cc    4-dimensional interpolation in the lookup table
-c
-         do ij = 1,nd
-           j(l,ij) =
-     $     poids(1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,ij)
-     $   + poids(1,1,1,2)*jphot(indtemp,indsza,indcol,indozo+1,ij)
-     $   + poids(1,1,2,1)*jphot(indtemp,indsza,indcol+1,indozo,ij)
-     $   + poids(1,1,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,ij)
-     $   + poids(1,2,1,1)*jphot(indtemp,indsza+1,indcol,indozo,ij)
-     $   + poids(1,2,1,2)*jphot(indtemp,indsza+1,indcol,indozo+1,ij)
-     $   + poids(1,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo,ij)
-     $   + poids(1,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,ij)
-     $   + poids(2,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,ij)
-     $   + poids(2,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo+1,ij)
-     $   + poids(2,1,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo,ij)
-     $   + poids(2,1,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,ij)
-     $   + poids(2,2,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,ij)
-     $   + poids(2,2,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,ij)
-     $   + poids(2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,ij)
-     $   + poids(2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,ij)
-         end do
-c
-cc    correction for sun distance
-c
-         do ij = 1,nd
-            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
-         end do 
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     end of loop over vertical levels
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      end do
-c
-      else
-c
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     night
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do ij = 1,nd
-         do l = 1,lswitch-1
-            j(l,ij) = 0.
-         end do
-      end do
-c
-      end if
-c
-      return
-      end
-c
-c*****************************************************************
-c
-      subroutine gcmtochim(zycol, lswitch, nesp, rm)
-c
-c*****************************************************************
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "callkeys.h"
-#include "tracer.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real zycol(nlayermx,nqmx)! species volume mixing ratio in the gcm
-c
-      integer nesp            ! number of species in the chemistry
-      integer lswitch         ! interface level between chemistries
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     outputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real rm(nlayermx,nesp)   ! species volume mixing ratio 
-c      
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer l,iq
-      logical,save :: firstcall=.true.
-      
-! Tracer indexes in the GCM:
-      integer,save :: g_co2=0
-      integer,save :: g_co=0
-      integer,save :: g_o=0
-      integer,save :: g_o1d=0
-      integer,save :: g_o2=0
-      integer,save :: g_o3=0
-      integer,save :: g_h=0
-      integer,save :: g_h2=0
-      integer,save :: g_oh=0
-      integer,save :: g_ho2=0
-      integer,save :: g_h2o2=0
-      integer,save :: g_n2=0
-      integer,save :: g_h2o_vap=0
-! 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_n2=13
-      integer,parameter :: i_hox=14
-      integer,parameter :: i_ox=15
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     tracers numbering in the gcm
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     co2      =  nqchem_min           
-c     co       =  nqchem_min + 1          
-c     o        =  nqchem_min + 2
-c     o(1d)    =  nqchem_min + 3
-c     o2       =  nqchem_min + 4
-c     o3       =  nqchem_min + 5
-c     h        =  nqchem_min + 6
-c     h2       =  nqchem_min + 7
-c     oh       =  nqchem_min + 8
-c     ho2      =  nqchem_min + 9
-c     h2o2     =  nqchem_min + 10
-c     n2       =  nqchem_min + 11
-c     ar       =  nqchem_min + 12
-c     h2o      =  nqmx
-c
-
-! Fist call initializations
-      if (firstcall) then
-        ! get the indexes of the tracers we'll need
-        do iq=1,nqmx
-          if (noms(iq).eq."co2") then
-            g_co2=iq
-          endif
-          if (noms(iq).eq."co") then
-            g_co=iq
-          endif
-          if (noms(iq).eq."o") then
-            g_o=iq
-          endif
-          if (noms(iq).eq."o1d") then
-            g_o1d=iq
-          endif
-          if (noms(iq).eq."o2") then
-            g_o2=iq
-          endif
-          if (noms(iq).eq."o3") then
-            g_o3=iq
-          endif
-          if (noms(iq).eq."h") then
-            g_h=iq
-          endif
-          if (noms(iq).eq."h2") then
-            g_h2=iq
-          endif
-          if (noms(iq).eq."oh") then
-            g_oh=iq
-          endif
-          if (noms(iq).eq."ho2") then
-            g_ho2=iq
-          endif
-          if (noms(iq).eq."h2o2") then
-            g_h2o2=iq
-          endif
-          if (noms(iq).eq."n2") then
-            g_n2=iq
-          endif
-          if (noms(iq).eq."h2o_vap") then
-            g_h2o_vap=iq
-          endif
-        enddo
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialise chemical species
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-       rm(l,i_co2)  = max(zycol(l, g_co2),     1.e-30)
-       rm(l,i_co)   = max(zycol(l, g_co),      1.e-30)
-       rm(l,i_o)    = max(zycol(l, g_o),       1.e-30)
-       rm(l,i_o1d)  = max(zycol(l, g_o1d),     1.e-30)
-       rm(l,i_o2)   = max(zycol(l, g_o2),      1.e-30)
-       rm(l,i_o3)   = max(zycol(l, g_o3),      1.e-30)
-       rm(l,i_h)    = max(zycol(l, g_h),       1.e-30)
-       rm(l,i_h2)   = max(zycol(l, g_h2),      1.e-30)
-       rm(l,i_oh)   = max(zycol(l, g_oh),      1.e-30)
-       rm(l,i_ho2)  = max(zycol(l, g_ho2),     1.e-30)
-       rm(l,i_h2o2) = max(zycol(l, g_h2o2),    1.e-30)
-       rm(l,i_n2)   = max(zycol(l, g_n2),      1.e-30)
-       rm(l,i_h2o)  = max(zycol(l, g_h2o_vap), 1.e-30)
-c 
-      end do 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     initialise chemical families                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      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
-c
-      return
-      end
-c
-c*****************************************************************
-c
-      subroutine chimtogcm(zycol, lswitch, nesp, rm)
-c
-c*****************************************************************
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "callkeys.h"
-#include"tracer.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer nesp               ! number of species in the chemistry
-      integer lswitch            ! interface level between chemistries
-c      
-      real rm(nlayermx,nesp)     ! species volume mixing ratio 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real zycol(nlayermx,nqmx)  ! species volume mixing ratio in the gcm
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer l,iq
-      logical,save :: firstcall=.true.
-      
-! Tracer indexes in the GCM:
-      integer,save :: g_co2=0
-      integer,save :: g_co=0
-      integer,save :: g_o=0
-      integer,save :: g_o1d=0
-      integer,save :: g_o2=0
-      integer,save :: g_o3=0
-      integer,save :: g_h=0
-      integer,save :: g_h2=0
-      integer,save :: g_oh=0
-      integer,save :: g_ho2=0
-      integer,save :: g_h2o2=0
-      integer,save :: g_n2=0
-      integer,save :: g_h2o_vap=0
-! 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_n2=13
-      integer,parameter :: i_hox=14
-      integer,parameter :: i_ox=15
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     tracers numbering in the gcm
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c     co2      =  nqchem_min
-c     co       =  nqchem_min + 1
-c     o        =  nqchem_min + 2
-c     o(1d)    =  nqchem_min + 3
-c     o2       =  nqchem_min + 4
-c     o3       =  nqchem_min + 5
-c     h        =  nqchem_min + 6
-c     h2       =  nqchem_min + 7
-c     oh       =  nqchem_min + 8
-c     ho2      =  nqchem_min + 9
-c     h2o2     =  nqchem_min + 10
-c     n2       =  nqchem_min + 11
-c     ar       =  nqchem_min + 12
-c     h2o      =  nqmx
-c      
-! Fist call initializations
-      if (firstcall) then
-        ! get the indexes of the tracers we'll need
-        do iq=1,nqmx
-          if (noms(iq).eq."co2") then
-            g_co2=iq
-          endif
-          if (noms(iq).eq."co") then
-            g_co=iq
-          endif
-          if (noms(iq).eq."o") then
-            g_o=iq
-          endif
-          if (noms(iq).eq."o1d") then
-            g_o1d=iq
-          endif
-          if (noms(iq).eq."o2") then
-            g_o2=iq
-          endif
-          if (noms(iq).eq."o3") then
-            g_o3=iq
-          endif
-          if (noms(iq).eq."h") then
-            g_h=iq
-          endif
-          if (noms(iq).eq."h2") then
-            g_h2=iq
-          endif
-          if (noms(iq).eq."oh") then
-            g_oh=iq
-          endif
-          if (noms(iq).eq."ho2") then
-            g_ho2=iq
-          endif
-          if (noms(iq).eq."h2o2") then
-            g_h2o2=iq
-          endif
-          if (noms(iq).eq."n2") then
-            g_n2=iq
-          endif
-          if (noms(iq).eq."h2o_vap") then
-            g_h2o_vap=iq
-          endif
-        enddo
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     save mixing ratios for the gcm
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-         zycol(l, g_co2)     = rm(l,i_co2) 
-         zycol(l, g_co)      = rm(l,i_co) 
-         zycol(l, g_o)       = rm(l,i_o) 
-         zycol(l, g_o1d)     = rm(l,i_o1d)
-         zycol(l, g_o2)      = rm(l,i_o2) 
-         zycol(l, g_o3)      = rm(l,i_o3) 
-         zycol(l, g_h)       = rm(l,i_h)  
-         zycol(l, g_h2)      = rm(l,i_h2) 
-         zycol(l, g_oh)      = rm(l,i_oh) 
-         zycol(l, g_ho2)     = rm(l,i_ho2) 
-         zycol(l, g_h2o2)    = rm(l,i_h2o2)
-         zycol(l, g_n2)      = rm(l,i_n2)
-         zycol(l, g_h2o_vap) = rm(l,i_h2o)
-      end do 
-c
-c  water ice: zycol(nqmx-1) has not been changed if iceparty on
-! not an issue anymore since we now use tracer names ...
-! and not relevant anymore as we removed iceparty: water=>iceparty
-c  
-      return
-      end
-c
-c*****************************************************************
-c
-      subroutine chemrates(lswitch, dens, press, t, icesurf,
-     $                     a001, a002, a003,
-     $                     b001, b002, b003, b004, b005, b006,
-     $                     c001, c002, c003, c004, c005, c006,
-     $                     c007, c008, c009, c010, c011, c012,
-     $                     c013, c014, c015, c016, c017, c018,
-     $                     d001, d002, d003,
-     $                     e001, e002,
-     $                     h001, h002)
-c
-c*****************************************************************
-c
-      implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:                                                    c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer lswitch                  ! interface level between chemistries
-
-      real dens(nlayermx)              ! density (cm-3)
-      real press(nlayermx)             ! pressure (hpa)
-      real t(nlayermx)                 ! temperature (k)
-      real icesurf(nlayermx)           ! ice surface area (cm^2/cm^3)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     outputs:                                                   c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
-      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
-     $     b004(nlayermx), b005(nlayermx), b006(nlayermx)
-      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
-     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
-     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
-     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
-     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
-     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
-      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
-      real e001(nlayermx), e002(nlayermx)
-      real h001(nlayermx), h002(nlayermx)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:                                                     c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real ak0, ak1, rate, xpo
-c
-      integer l
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     compute reaction rates from JPL 1997, otherwise mentioned
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      do l = 1,lswitch-1
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        oxygen compounds
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      a001: o + o2 + co2 -> o3 + co2
-c
-c        jpl 2003
-c
-         a001(l) = 2.5
-     $             *6.0e-34*(t(l)/300.)**(-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
-         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)
-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
-         b001(l) = 7.4e-11*exp(120./t(l))
-c
-ccc      b002: o(1d) + h2o  -> oh + oh
-c
-c        jpl 2003
-c
-         b002(l) = 2.2e-10
-c    
-ccc      b003: o(1d) + h2  -> oh + h
-c
-c        jpl 2003
-c
-         b003(l) = 1.1e-10
-c
-c        nair et al., 1994
-c
-c        b003(l) = 1.0e-10
-c    
-ccc      b004: o(1d) + o2  -> o + o2
-c
-c        jpl 2003
-c
-         b004(l) = 3.2e-11*exp(70./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
-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        c001(l) = 0.75*3.0e-11*exp(200./t(l))
-c
-ccc      c002: o + oh -> o2 + h
-c
-c        jpl 2003
-c 
-         c002(l) = 2.2e-11*exp(120./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 
-         c004(l) = 8.1e-11*0.90
-c
-ccc      c005: h + ho2 -> h2 + o2
-c
-c        jpl 2003
-c 
-         c005(l) = 8.1e-11*0.08
-c
-ccc      c006: h + ho2 -> h2o + o
-c
-c        jpl 2003
-c 
-         c006(l) = 8.1e-11*0.02
-c
-ccc      c007: oh + ho2 -> h2o + o2
-c
-c        jpl 2003
-c 
-         c007(l) = 4.8e-11*exp(250./t(l))
-c        c007(l) = 0.75*4.8e-11*exp(250./t(l))
-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 
-         c009(l) = 2.9e-12*exp(-160./t(l))
-c
-ccc      c010: oh + h2 -> h2o + h
-c
-c        jpl 2003
-c 
-         c010(l) = 5.5e-12*exp(-2000./t(l))
-c
-ccc      c011: h + o2 + co2 -> ho2 + co2
-c
-c        jpl 2003
-c 
-         ak0 = 2.5*5.7e-32*(t(l)/300.)**(-1.6)
-         ak1 = 7.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)
-         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 
-         c013(l) = 4.2e-12*exp(-240./t(l))
-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 2003
-c
-         c016(l) = 2.5*1.7e-33
-     $              *exp(1000./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
-         c018(l) = 2.5*8.85e-33*(t(l)/298.)**(-0.6)*dens(l)
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        nitrogen compounds
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      d001: no2 + o -> no + o2
-c
-c        jpl 2003
-c
-         d001(l) = 5.6e-12*exp(180./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 2003
-c
-         d003(l) = 3.5e-12*exp(250./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
-         e001(l) = 1.57e-13 + 3.54e-33*dens(l)
-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
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c        heterogenous chemistry
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-ccc      h001: ho2 + ice surface -> (ho2)ads
-ccc           (ho2)ads + oh -> h2o + o2
-ccc      total : oh + ho2 -> h2o + o2, first order
-ccc         K = h001*[ho2],
-ccc      h001 = (S*v*gamma1)/4. (s-1), v=100*sqrt(3kT/m_ho2) (cm s-1)
-c
-c        cooper and abbatt, 1996
-c
-         h001(l) = icesurf(l) * 25.*sqrt(3*8.31*t(l)/33.) * 0.025
-c         h001(l) = 0.0e0
-c
-ccc      h002: oh + ice surface -> (oh)ads
-ccc           (oh)ads + ho2 -> h2o + o2
-ccc      total : oh + ho2 -> h2o + o2, first order
-ccc         K = h002*[oh],
-ccc      h002 = (S*v*gamma2)/4. (s-1), v=100*sqrt(3kT/m_oh) (cm s-1)
-c
-c        cooper and abbatt, 1996
-c
-         h002(l) = icesurf(l) * 25.*sqrt(3*8.31*t(l)/17.) * 0.03
-c         h002(l) = 0.0e0
-c
-c        write(*,*) "level ",l," / icesurf=",icesurf(l),
-c     $       " / 0.25vg1=",(25.*sqrt(3*8.31*t(l)/33.) * 0.025),
-c     $       " / 0.25vg2=",(25.*sqrt(3*8.31*t(l)/17.) * 0.03)
-c
-      end do
-c
-      return
-      end
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 334)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 334)
@@ -0,0 +1,2007 @@
+c*****************************************************************
+c
+c     Photochemical routine 
+c
+c     Author: Franck Lefevre
+c     ------
+c
+c     Version: 17/03/2011
+c
+c*****************************************************************
+c
+      subroutine photochemistry(lswitch, zycol, sza, ptimestep, press,
+     $                          temp, dens, dist_sol, surfdust1d,
+     $                          surfice1d, jo3, tau)
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "chimiedata.h" 
+#include "callkeys.h"
+#include "tracer.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     input/output:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      real zycol(nlayermx,nqmx)  ! chemical species volume mixing ratio
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      real sza                   ! solar zenith angle (deg)
+      real ptimestep             ! physics timestep (s)
+      real press(nlayermx)       ! pressure (hpa)
+      real temp(nlayermx)        ! temperature (k)
+      real dens(nlayermx)        ! density (cm-3)
+      real dist_sol              ! sun distance (au) 
+      real surfdust1d(nlayermx)  ! dust surface area (cm2/cm3)
+      real surfice1d(nlayermx)   ! ice surface area (cm2/cm3)
+      real tau                   ! optical depth at 7 hpa
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     output:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      real jo3(nlayermx)  ! photodissociation rate o3 -> o1d
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer phychemrat         ! ratio physics timestep/chemistry timestep
+      integer istep
+      integer i_co2,i_o3,j_o3_o1d,lev
+      integer nesp
+      integer lswitch
+c
+      parameter (nesp = 16)      ! number of species in the chemistry
+c
+      real ctimestep             ! chemistry timestep (s) 
+      real rm(nlayermx,nesp)     ! species volume mixing ratio 
+      real j(nlayermx,nd)        ! interpolated photolysis rates (s-1)
+      real rmco2(nlayermx)       ! co2 mixing ratio
+      real rmo3(nlayermx)        ! ozone mixing ratio
+c
+c     reaction rates
+c
+      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
+      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
+     $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
+     $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
+      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
+     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
+     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
+     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
+     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
+     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
+      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
+      real e001(nlayermx), e002(nlayermx), e003(nlayermx)
+      real h001(nlayermx), h002(nlayermx), h003(nlayermx),
+     $     h004(nlayermx), h005(nlayermx)
+      real t001(nlayermx), t002(nlayermx), t003(nlayermx)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     ctimestep : chemistry timestep (s)                         c
+c                 10 minutes                                     c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      ctimestep = 10.*60.
+c
+      phychemrat = int(ptimestep/ctimestep)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     initialisation of chemical species and families            c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      call gcmtochim(zycol, lswitch, nesp, rm) 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     compute reaction rates                                     c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      call chemrates(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)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     interpolation of photolysis rates in the lookup table      c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      i_co2 = 1
+      i_o3  = 6 
+c
+      do lev = 1,lswitch-1
+         rmco2(lev) = rm(lev,i_co2)
+         rmo3(lev)  = rm(lev, i_o3)
+      end do
+c
+      call phot(lswitch, press, temp, sza, tau, dist_sol, 
+     $          rmco2, rmo3, j)
+c
+      j_o3_o1d = 5
+c
+      do lev = 1,lswitch-1
+         jo3(lev) = j(lev,j_o3_o1d)
+      end do
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     loop over time within the physical timestep                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do istep = 1,phychemrat
+         call chimie(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
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     save chemical species for the gcm                          c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      call chimtogcm(zycol, lswitch, nesp, rm) 
+c
+      return
+      end
+c
+c*****************************************************************
+c
+      subroutine chimie(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
+c*****************************************************************
+c
+      implicit none 
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "chimiedata.h" 
+#include "callkeys.h" 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     input/output:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      integer lswitch         ! interface level between chemistries
+      integer nesp            ! number of species
+c
+      real rm(nlayermx,nesp)  ! volume mixing ratios 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real dens(nlayermx)     ! density (cm-3)
+      real dt                 ! chemistry timestep (s) 
+      real j(nlayermx,nd)     ! interpolated photolysis rates (s-1) 
+      real press(nlayermx)    ! pressure (hpa)
+      real t(nlayermx)        ! temperature (k)
+      real sza                ! solar zenith angle (deg)
+      real dist_sol           ! sun distance (au) 
+c
+c     reaction rates
+c
+      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
+      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
+     $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
+     $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
+      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
+     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
+     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
+     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
+     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
+     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
+      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
+      real e001(nlayermx), e002(nlayermx), e003(nlayermx)
+      real h001(nlayermx), h002(nlayermx), h003(nlayermx),
+     $     h004(nlayermx), h005(nlayermx)
+      real t001(nlayermx), t002(nlayermx), t003(nlayermx)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real hetero_ice, hetero_dust
+c
+      integer iesp, iter, l, niter
+      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 
+c
+      parameter (hetero_ice  = 1.)   ! switch for het. chem. on ice clouds
+      parameter (hetero_dust = 0.)   ! switch for het. chem. on dust
+                                     ! hetero_dust = 0. advised for the moment
+c
+      parameter (niter = 5)          ! iterations in the chemical scheme
+c
+      real cc0(nlayermx,nesp)        ! initial number density (cm-3)
+      real cc(nlayermx,nesp)         ! final number density (cm-3)
+      real nox(nlayermx)             ! nox number density (cm-3)
+      real no(nlayermx)              ! no  number density (cm-3)
+      real no2(nlayermx)             ! no2 number density (cm-3)
+      real production(nlayermx,nesp) ! production rate
+      real loss(nlayermx,nesp)       ! loss rate
+c
+      real ro_o3, rh_ho2, roh_ho2, rno2_no
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     tracer numbering in the chemistry
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     numbering of photolysis rates
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     volume mixing ratio -> density conversion
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     co2 and nox number densities (cm-3)    
+c      
+c     nox mixing ratio: 6.e-10 (yung and demore, 1999)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+         nox(l) = 6.e-10*dens(l)
+      end do 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     loop over iterations in the chemical scheme
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do iter = 1,niter
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        nox species
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        no2/no
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+c
+         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)))
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        no
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         no(l) = nox(l)/(1. + rno2_no)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        no2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         no2(l) = no(l)*rno2_no
+c
+      end do
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        hox species
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        photochemical equilibrium for oh and ho2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        h/ho2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+c
+         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
+     $           )
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        oh/ho2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        h
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         cc(l,i_h) = cc(l,i_hox)
+     $                 /(1. + (1. + roh_ho2)/rh_ho2)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        ho2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         cc(l,i_ho2) = cc(l,i_h)/rh_ho2
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        oh
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
+c
+      end do
+c
+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
+c        
+      if (sza .le. 95.) then
+c
+         do l = 1,lswitch-1
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o(1d)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            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))
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o/o3
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            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))
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o3
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            cc(l,i_o) = cc(l,i_o3)*ro_o3
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        ox = o + o3
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            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)
+c
+            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)
+c
+            loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
+c
+         end do
+c
+      else
+c
+         do l = 1,lswitch-1
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o(1d)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            cc(l,i_o1d) = 0.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o3
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            production(l,i_o3) = a001(l)*cc(l,i_o2)*cc(l,i_o)
+c
+            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)
+c
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+            production(l,i_o) = c006(l)*cc(l,i_h)*cc(l,i_ho2)
+     $                        + c013(l)*cc(l,i_oh)*cc(l,i_oh)
+c
+            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)
+c
+         end do
+      end if
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     other species
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        co2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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
+c
+         loss(l,i_co2) = j(l,j_co2_o)
+     $                 + j(l,j_co2_o1d)
+     $                 + t003(l)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        co
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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)
+c
+         loss(l,i_co) = e001(l)*cc(l,i_oh)
+     $                + e002(l)*cc(l,i_o)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        o2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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)
+c
+         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)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        h2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         production(l,i_h2) = c005(l)*cc(l,i_h)*cc(l,i_ho2)
+     $                      + c018(l)*cc(l,i_h)*cc(l,i_h)
+c
+         loss(l,i_h2) = b003(l)*cc(l,i_o1d)
+     $                + c010(l)*cc(l,i_oh)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        h2o
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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
+c
+         loss(l,i_h2o) = j(l,j_h2o)
+     $                 + b002(l)*cc(l,i_o1d)
+     $                 + t001(l)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        h2o2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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
+c
+         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
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        hox = h + oh + ho2
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         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)
+c
+         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
+c
+         loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        ch4
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+         production(l,i_ch4) = 0.
+c
+         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)
+c
+      end do
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     update number densities
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     long-lived species
+c
+      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
+c     ox species
+c
+      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 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     end of loop over iterations
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      end do
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     density -> volume mixing ratio conversion
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c
+      return
+      end
+c
+c*****************************************************************
+c
+      subroutine phot(lswitch, press, temp, sza, tauref, dist_sol, 
+     $                rmco2, rmo3, j)
+c
+c*****************************************************************
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "chimiedata.h" 
+#include "comcstfi.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c       
+      integer lswitch         ! interface level between chemistries
+      real press(nlayermx)             ! pressure (hPa)
+      real temp(nlayermx)              ! temperature (K)
+      real sza                         ! solar zenith angle (deg)
+      real tauref                      ! optical depth at 7 hpa
+      real dist_sol                    ! sun distance (AU) 
+      real rmco2(nlayermx)             ! co2 mixing ratio
+      real rmo3(nlayermx)              ! ozone mixing ratio
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     output:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real j(nlayermx,nd)              ! interpolated photolysis rates (s-1) 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer icol, ij, indsza, indtau, indcol, indozo, indtemp, 
+     $        iozo, isza, itau, it, l
+c
+      real col(nlayermx)               ! overhead air column   (molecule cm-2)
+      real colo3(nlayermx)             ! 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(nlayermx)
+      real tau
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     day/night criterion
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      if (sza .le. 95.) then 
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     temperatures at 1.9 hPa in lookup table
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      table_temp(1) = 226.2
+      table_temp(2) = 206.2
+      table_temp(3) = 186.2
+      table_temp(4) = 169.8
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     interpolation in solar zenith angle
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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))
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     interpolation in dust (tau)
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      tau = min(tauref, tautab(ntau))
+      tau = max(tau, tautab(1))
+c
+      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))
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     co2 and ozone columns
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     co2 column at model top (molecule.cm-2)
+c
+      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.
+     $                 /(mugaz*g)
+c
+c     ozone column at model top
+c
+      colo3(lswitch-1) = 0.
+c
+c     co2 and ozone columns for other levels (molecule.cm-2)
+c
+      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
+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.
+c
+      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
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     temperature dependence
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     1) search for temperature at 1.9 hPa (tref): vertical interpolation
+c
+      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
+c     2) interpolation in temperature
+c
+      tref = min(tref,table_temp(1))
+      tref = max(tref,table_temp(ntemp))
+c
+      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
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     loop over vertical levels
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+c
+c     interpolation in air column
+c
+         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
+c
+cc    interpolation in ozone column
+c
+         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.)
+c
+cc    4-dimensional interpolation weights
+c
+c     poids(temp,sza,co2,o3,tau)
+c
+         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)
+c
+         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
+c
+cc    4-dimensional interpolation in the lookup table
+c
+         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)
+c
+     $   + 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
+c
+cc    correction for sun distance
+c
+         do ij = 1,nd
+            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
+         end do 
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     end of loop over vertical levels
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      end do
+c
+      else
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     night
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do ij = 1,nd
+         do l = 1,lswitch-1
+            j(l,ij) = 0.
+         end do
+      end do
+c
+      end if
+c
+      return
+      end
+c
+c*****************************************************************
+c
+      subroutine gcmtochim(zycol, lswitch, nesp, rm)
+c
+c*****************************************************************
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "callkeys.h"
+#include "tracer.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      real zycol(nlayermx,nqmx)! species volume mixing ratio in the gcm
+c
+      integer nesp             ! number of species in the chemistry
+      integer lswitch          ! interface level between chemistries
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     outputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real rm(nlayermx,nesp)   ! species volume mixing ratio 
+c      
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer l,iq
+      logical,save :: firstcall = .true.
+      
+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
+c      
+c     first call initializations
+c
+      if (firstcall) then
+
+c       identify the indexes of the tracers we need
+
+        if (igcm_co2.eq.0) then
+          write(*,*) "concentrations: Error; no CO2 tracer !!!"
+          stop
+        endif
+        if (igcm_co.eq.0) then
+          write(*,*) "concentrations: Error; no CO tracer !!!"
+          stop
+        endif
+        if (igcm_o.eq.0) then
+          write(*,*) "concentrations: Error; no O tracer !!!"
+          stop
+        endif
+        if (igcm_o1d.eq.0) then
+          write(*,*) "concentrations: Error; no O1D tracer !!!"
+          stop
+        endif
+        if (igcm_o2.eq.0) then
+          write(*,*) "concentrations: Error; no O2 tracer !!!"
+          stop
+        endif
+        if (igcm_o3.eq.0) then
+          write(*,*) "concentrations: Error; no O3 tracer !!!"
+          stop
+        endif
+        if (igcm_h.eq.0) then
+          write(*,*) "concentrations: Error; no H tracer !!!"
+          stop
+        endif
+        if (igcm_h2.eq.0) then
+          write(*,*) "concentrations: Error; no H2 tracer !!!"
+          stop
+        endif
+        if (igcm_oh.eq.0) then
+          write(*,*) "concentrations: Error; no OH tracer !!!"
+          stop
+        endif
+        if (igcm_ho2.eq.0) then
+          write(*,*) "concentrations: Error; no HO2 tracer !!!"
+          stop
+        endif
+        if (igcm_h2o2.eq.0) then
+          write(*,*) "concentrations: Error; no H2O2 tracer !!!"
+          stop
+        endif
+        if (igcm_ch4.eq.0) then
+          write(*,*) "concentrations: Error; no CH4 tracer !!!"
+          stop
+        endif
+        if (igcm_n2.eq.0) then
+          write(*,*) "concentrations: Error; no N2 tracer !!!"
+          stop
+        endif
+        if (igcm_ar.eq.0) then
+          write(*,*) "concentrations: Error; no AR tracer !!!"
+          stop
+        endif
+        if (igcm_h2o_vap.eq.0) then
+          write(*,*) "concentrations: Error; no water vapor tracer !!!"
+          stop
+        endif
+        firstcall = .false.
+      end if ! of if (firstcall)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     initialise chemical species
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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_ch4)  = max(zycol(l, igcm_ch4),     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 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     initialise chemical families                     c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c
+      return
+      end
+c
+c*****************************************************************
+c
+      subroutine chimtogcm(zycol, lswitch, nesp, rm)
+c
+c*****************************************************************
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "callkeys.h"
+#include "tracer.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer nesp               ! number of species in the chemistry
+      integer lswitch            ! interface level between chemistries
+c      
+      real rm(nlayermx,nesp)     ! species volume mixing ratio 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     output:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      
+      real zycol(nlayermx,nqmx)  ! species volume mixing ratio in the gcm
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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
+c      
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     save mixing ratios for the gcm
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      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_ch4)     = rm(l,i_ch4)
+         zycol(l, igcm_n2)      = rm(l,i_n2)
+         zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
+      end do 
+c
+      return
+      end
+c
+c*****************************************************************
+c
+      subroutine chemrates(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
+c*****************************************************************
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:                                                    c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer lswitch                  ! interface level between chemistries
+
+      real dens(nlayermx)              ! density (cm-3)
+      real press(nlayermx)             ! pressure (hpa)
+      real t(nlayermx)                 ! temperature (k)
+      real surfdust1d(nlayermx)        ! dust surface area (cm^2/cm^3)
+      real surfice1d(nlayermx)         ! ice surface area (cm^2/cm^3)
+      real tribo                       ! switch for triboelectricity
+      real tau                         ! dust opacity at 7 hpa
+c
+      parameter (tribo   = 0.)         ! switch for triboelectricity
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     outputs:                                                   c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real a001(nlayermx), a002(nlayermx), a003(nlayermx)
+      real b001(nlayermx), b002(nlayermx), b003(nlayermx), 
+     $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
+     $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
+      real c001(nlayermx), c002(nlayermx), c003(nlayermx), 
+     $     c004(nlayermx), c005(nlayermx), c006(nlayermx), 
+     $     c007(nlayermx), c008(nlayermx), c009(nlayermx), 
+     $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
+     $     c013(nlayermx), c014(nlayermx), c015(nlayermx), 
+     $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
+      real d001(nlayermx), d002(nlayermx), d003(nlayermx)
+      real e001(nlayermx), e002(nlayermx), e003(nlayermx)
+      real h001(nlayermx), h002(nlayermx), h003(nlayermx),
+     $     h004(nlayermx), h005(nlayermx)
+      real t001(nlayermx), t002(nlayermx), t003(nlayermx)
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:                                                     c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      real ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
+      real ef, efmax, lossh2o, lossch4, lossco2
+c
+      integer l
+      real k1a, k1b, k1a0, k1b0, k1ainf
+      real x, y, fc, fx
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     compute reaction rates
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      do l = 1,lswitch-1
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        oxygen compounds
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+ccc      a001: o + o2 + co2 -> o3 + co2
+c
+c        jpl 2003
+c
+         a001(l) = 2.5
+     $             *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 2003
+c
+         b003(l) = 1.1e-10
+c
+c        nair et al., 1994
+c
+c        b003(l) = 1.0e-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 2003
+c 
+         c002(l) = 2.2e-11*exp(120./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 2003
+c 
+c        ak0 = 2.5*5.7e-32*(t(l)/300.)**(-1.6)
+c        ak1 = 7.5e-11*(t(l)/300.)**(0.0)
+c
+c        jpl 2006
+c
+         ak0 = 2.5*4.4e-32*(t(l)/300.)**(-1.3)
+         ak1 = 4.7e-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 2003
+c
+         c016(l) = 2.5*1.7e-33
+     $              *exp(1000./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 2003
+c
+         d003(l) = 3.5e-12*exp(250./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
Index: trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F	(revision 334)
+++ trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F	(revision 334)
@@ -0,0 +1,79 @@
+      subroutine read_phototable
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c   subject:
+c   --------
+c
+c   read photolysis lookup table
+c
+c   VERSION: 17/03/2011
+c
+c   Author:   Franck Lefevre
+c
+c   Arguments:
+c   ----------
+c
+c   The output variable is jphot and is put in common chimiedata.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "chimiedata.h"
+#include "datafile.h"
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     local:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      integer fic, ij, iozo, isza, itemp, iz, itau, ierr
+c
+      real xsza
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      fic = 81
+
+      open(fic, form = 'formatted', status = 'old',
+     &file =trim(datafile)//'/jmars.20101220',iostat=ierr)
+
+      IF (ierr.NE.0) THEN
+        write(*,*)'Error : cannot open photolysis lookup table'
+        write(*,*)'(in aeronomars/init_chemistry.F)'
+        write(*,*)'It should be in :',trim(datafile),'/'
+        write(*,*)'1) You can change this directory address in '
+        write(*,*)'   file phymars/datafile.h'
+        STOP
+      ENDIF
+
+c
+      print*, 'read photolysis lookup table...'
+c
+      do itau = 1,ntau
+         do itemp = 1,ntemp
+            do iozo = 1,nozo
+               do isza = 1,nsza
+                  do iz = 200,0,-1
+                     read(fic,*) colairtab(iz), xsza, table_ozo(iozo)
+c                    print*, colairtab(iz), xsza, table_ozo(iozo)
+                     read(fic,'(7e11.4)')
+     $               (jphot(itemp,isza,iz,iozo,itau,ij), ij= 1,nd)
+                     do ij = 1,nd
+                        if (jphot(itemp,isza,iz,iozo,itau,ij)
+     $                     .eq. 1.e-30) then
+                           jphot(itemp,isza,iz,iozo,itau,ij) = 0.
+                        end if
+                     end do
+                  end do
+               end do
+            end do
+         end do
+      end do
+c
+      print*, 'lookup table...ok'
+      close(fic)
+c
+      return
+      end
Index: trunk/LMDZ.MARS/libf/dyn3d/inidissip.F
===================================================================
--- trunk/LMDZ.MARS/libf/dyn3d/inidissip.F	(revision 333)
+++ trunk/LMDZ.MARS/libf/dyn3d/inidissip.F	(revision 334)
@@ -62,4 +62,14 @@
              delta=30.! Intervalle (km) pour le changement mid /up
          end if
+
+         !!!!! reglages 35 niveaux FL
+         !fac_mid=3
+         !fac_up=30
+         !startalt=70.
+         !delta=20.
+         !!!! avec dans run.def
+         !!- mode_sponge=3
+         !!- idissip/tetagdiv/tetagrot/tetagtemp = 1/2000/3000/3000
+
          middle=startalt+delta/2
          write(*,*) 'Dissipation : '
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 333)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 334)
@@ -198,4 +198,6 @@
       REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
                                      !   of the size distribution
+      REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry)
+      REAL surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3, if photochemistry)
 
 c     Variables used by the slope model
@@ -306,4 +308,5 @@
       real rho(ngridmx,nlayermx)  ! density
       real vmr(ngridmx,nlayermx)  ! volume mixing ratio
+      !real colden(ngridmx,nqmx)   ! vertical column                      !FL
       REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
       REAL icetot(ngridmx)        ! Total mass of water ice (kg/m2)
@@ -918,5 +921,7 @@
      &                pq,pdq,zdqcloud,zdqscloud,zdtcloud,
      &                nq,naerkind,tau,
-     &                ccn,rdust,rice,nuice)
+     &                ccn,rdust,rice,nuice,
+     &                surfdust,surfice)
+
            if (activice) then
 c Temperature variation due to latent heat release
@@ -955,8 +960,8 @@
 c        --------------
          IF (photochem .or. thermochem) then
-          call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
-     &      zzlay,zday,pq,pdq,rice,
-     &      zdqchim,zdqschim,zdqcloud,zdqscloud)
-!NB: Photochemistry includes condensation of H2O2
+            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
+     $                   zdqcloud,zdqscloud,tauref,co2ice,
+     $                   pu,pdu,pv,pdv,surfdust,surfice)
 
            ! increment values of tracers:
@@ -1409,8 +1414,20 @@
      .                     "Volume mixing ratio","mol/mol",3,vmr)
                    endif
-                enddo
-             endif ! of if (thermochem.or.photochem)
-
-           endif ! of if (tracer)
+!                   do ig = 1,ngrid
+!                      colden(ig,iq) = 0.                             !FL
+!                   end do
+!                   do l=1,nlayer                                     !FL
+!                      do ig=1,ngrid                                  !FL
+!                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) !FL
+!     $                                  *(pplev(ig,l)-pplev(ig,l+1)) !FL
+!     $                                  *6.022e22/(mmol(iq)*g)       !FL
+!                      end do                                         !FL
+!                   end do                                            !FL
+!                   call wstats(ngrid,"c_"//trim(noms(iq)),           !FL
+!     $                         "column","mol cm-2",2,colden(1,iq))   !FL
+                end do
+             end if ! of if (thermochem.or.photochem)
+
+           end if ! of if (tracer)
 
            IF(lastcall) THEN
@@ -1533,4 +1550,13 @@
          call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
      &                  co2col)
+!!!!! FL
+!            do iq = 1,nq
+!               if (noms(iq) .ne. "dust_mass" .and.
+!     $             noms(iq) .ne. "dust_number") then
+!               call writediagfi(ngrid,"c_"//trim(noms(iq)),         
+!     $                         "column","mol cm-2",2,colden(1,iq))
+!               end if
+!            end do
+!!!!! FL
          endif ! of if (tracer.and.(igcm_co2.ne.0))
 
Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 333)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 334)
@@ -3,5 +3,6 @@
      &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
      &                nq,naersize,tau,
-     &                ccn,rdust,rice,nuice)
+     &                ccn,rdust,rice,nuice,
+     &                surfdust,surfice)
       IMPLICIT NONE
 
@@ -70,4 +71,7 @@
                                !   of the size distribution
 
+      REAL surfdust(ngrid,nlay)   ! dust surface area (micron2/cm3, if photochemistry)
+      REAL surfice(ngrid,nlay)    !  ice surface area (micron2/cm3, if photochemistry)
+
 c   local:
 c   ------
@@ -85,4 +89,5 @@
       REAL masse (ngridmx,nlayermx)
       REAL epaisseur (ngridmx,nlayermx)
+      REAL npart(ngridmx,nlayermx)    !  Cloud condensation nuclei (#.cm-3)
 !      REAL rfinal        ! Ice crystal radius after condensation(m)
       REAL*8 seq           ! Equilibrium saturation ration (accounting for curvature effect)
@@ -179,4 +184,30 @@
       enddo ! of dol=1,nlay
 
+      if (photochem) then
+c        computation of dust and ice surface area (micron2/cm3)
+c        for heterogeneous chemistry
+
+         do l = 1,nlay
+            do ig = 1,ngrid
+c                                                                      
+c              npart: number density of ccn in #/cm3    
+c                                                     
+               npart(ig,l) = 1.e-6*ccn(ig,l)  
+     $                       *masse(ig,l)/epaisseur(ig,l) 
+c
+c              dust and ice surface area
+c                                                                  
+               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
+c                                                                  
+               if (rice(ig,l) .ge. rdust(ig,l)) then                   
+                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
+                  surfdust(ig,l) = 0.
+               else                                                    
+                  surfice(ig,l) = 0.                                 
+               end if                                             
+            end do      ! of do ig=1,ngrid
+         end do         ! of do l=1,nlay
+      end if            ! of photochem
+c
       pdqscloud(1:ngrid,1:nq)=0
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
