Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 1124)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 1125)
@@ -13,4 +13,5 @@
                             igcm_noplus, igcm_n2plus, igcm_hplus,         &
                             igcm_hco2plus, igcm_elec, mmol
+
       use conc_mod, only: mmean ! mean molecular mass of the atmosphere
 
@@ -33,4 +34,5 @@
 !    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
 !    update 16/03/2012 optimization (Franck Lefevre)
+!    update 11/12/2013 optimization (Franck Lefevre)
 !
 !   Arguments:
@@ -74,7 +76,7 @@
 !     input:
 
-      integer,intent(in) :: ngrid ! number of atmospheric columns
-      integer,intent(in) :: nlayer ! number of atmospheric layers
-      integer,intent(in) :: nq ! number of tracers
+      integer,intent(in) :: ngrid    ! number of atmospheric columns
+      integer,intent(in) :: nlayer   ! number of atmospheric layers
+      integer,intent(in) :: nq       ! number of tracers
       real :: ptimestep
       real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
@@ -88,11 +90,11 @@
       real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
       real :: pdv(ngrid,nlayer)      ! v component tendency
-      real :: dist_sol                   ! distance of the sun (AU)
-      real :: mu0(ngrid)               ! cos of solar zenith angle (=1 when sun at zenith)
-      real :: pq(ngrid,nlayer,nq)  ! tracers mass mixing ratio
-      real :: pdq(ngrid,nlayer,nq) ! previous tendencies
-      real :: zday                       ! date (time since Ls=0, in martian days)
-      real :: tauref(ngrid)            ! optical depth at 7 hPa
-      real :: co2ice(ngrid)            ! co2 ice surface layer (kg.m-2)
+      real :: dist_sol               ! distance of the sun (AU)
+      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
+      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
+      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
+      real :: zday                   ! date (time since Ls=0, in martian days)
+      real :: tauref(ngrid)          ! optical depth at 7 hPa
+      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
       real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
       real :: surfice(ngrid,nlayer)  !  ice surface area (m2/m3)
@@ -100,7 +102,7 @@
 !     output:
 
-      real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
+      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
       real :: dqschim(ngrid,nq)         ! tendencies on qsurf 
-      real :: dqcloud(ngrid,nlayer,nq)! tendencies on pq due to condensation
+      real :: dqcloud(ngrid,nlayer,nq)  ! tendencies on pq due to condensation
       real :: dqscloud(ngrid,nq)        ! tendencies on qsurf 
 
@@ -109,5 +111,5 @@
       integer,save :: nbq                   ! number of tracers used in the chemistry
       integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
-      integer :: iloc(1)            ! index of major species
+      integer :: iloc(1)                    ! index of major species
       integer :: ig,l,i,iq,iqmax
       integer :: foundswitch, lswitch
@@ -147,29 +149,29 @@
 
       real    :: latvl1, lonvl1
-      real    :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
-                                           ! new mole fraction after
+      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
+                                       ! new mole fraction after
       real    :: zt(ngrid,nlayer)      ! temperature
       real    :: zu(ngrid,nlayer)      ! u component of the wind
       real    :: zv(ngrid,nlayer)      ! v component of the wind
-      real    :: taucol                    ! optical depth at 7 hPa
+      real    :: taucol                ! optical depth at 7 hPa
 
       logical,save :: firstcall = .true.
-      logical,save :: depos = .false.      ! switch for dry deposition
+      logical,save :: depos = .false.  ! switch for dry deposition
 
 !     for each column of atmosphere:
 
-      real :: zpress(nlayer)       !  Pressure (mbar)
-      real :: zdens(nlayer)        !  Density  (cm-3)
-      real :: ztemp(nlayer)        !  Temperature (K)
-      real :: zlocal(nlayer)       !  Altitude (km)
-      real :: zycol(nlayer,nq)   !  Composition (mole fractions)
-      real :: szacol                 !  Solar zenith angle
-      real :: surfice1d(nlayer)    !  Ice surface area (cm2/cm3)
-      real :: surfdust1d(nlayer)   !  Dust surface area (cm2/cm3)
-      real :: jo3(nlayer)          !  Photodissociation rate O3->O1D (s-1)
+      real :: zpress(nlayer)       ! Pressure (mbar)
+      real :: zdens(nlayer)        ! Density  (cm-3)
+      real :: ztemp(nlayer)        ! Temperature (K)
+      real :: zlocal(nlayer)       ! Altitude (km)
+      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
+      real :: szacol               ! Solar zenith angle
+      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
+      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
+      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
 
 !     for output:
 
-      logical :: output                 ! to issue calls to writediagfi and stats
+      logical :: output             ! to issue calls to writediagfi and stats
       parameter (output = .true.)
       real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
@@ -645,5 +647,6 @@
 
          if (photochem) then
-            call photochemistry(lswitch,zycol,szacol,ptimestep,    &
+            call photochemistry(nlayer,nq,                         &
+                                lswitch,zycol,szacol,ptimestep,    &
                                 zpress,ztemp,zdens,dist_sol,       &
                                 surfdust1d,surfice1d,jo3,taucol)
Index: /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1124)
+++ /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1125)
@@ -6,84 +6,86 @@
 c     ------
 c
-c     Version: 17/03/2011
+c     Version: 11/12/2013
 c
 c*****************************************************************
-c
-      subroutine photochemistry(lswitch, zycol, sza, ptimestep, 
+
+      subroutine photochemistry(nlayer, nq,
+     $                  lswitch, zycol, sza, ptimestep, 
      $                  press,temp, dens, dist_sol, surfdust1d,
      $                  surfice1d, jo3, tau)
-c
-      use tracer_mod, only: nqmx
+
       implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "chimiedata.h" 
 #include "callkeys.h"
 !#include "tracer.h"
-c
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer, intent(in) :: nq     ! number of tracers
+      
+      real :: sza                   ! solar zenith angle (deg)
+      real :: ptimestep             ! physics timestep (s)
+      real :: press(nlayer)         ! pressure (hpa)
+      real :: temp(nlayer)          ! temperature (k)
+      real :: dens(nlayer)          ! density (cm-3)
+      real :: dist_sol              ! sun distance (au) 
+      real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
+      real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
+      real :: tau                   ! optical depth at 7 hpa
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     input/output:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-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
+      
+      real :: zycol(nlayer,nq)      ! chemical species volume mixing ratio
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     output:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real jo3(nlayermx)  ! photodissociation rate o3 -> o1d
-c
+      
+      real :: jo3(nlayer)  ! photodissociation rate o3 -> o1d
+
 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 stimestep             ! standard timestep for the chemistry (s)
-      real ctimestep             ! real timestep for the chemistry (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
+
+      integer, parameter :: nesp = 16 ! number of species in the chemistry
+
+      integer :: phychemrat       ! ratio physics timestep/chemistry timestep
+      integer :: istep
+      integer :: i_co2,i_o3,j_o3_o1d,lev
+      integer :: lswitch
+
+      real :: stimestep           ! standard timestep for the chemistry (s)
+      real :: ctimestep           ! real timestep for the chemistry (s) 
+      real :: rm(nlayer,nesp)     ! species volume mixing ratio 
+      real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
+      real :: rmco2(nlayer)       ! co2 mixing ratio
+      real :: rmo3(nlayer)        ! ozone mixing ratio
+
 c     reaction rates
-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
+
+      real :: a001(nlayer), a002(nlayer), a003(nlayer)
+      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
+     $        b004(nlayer), b005(nlayer), b006(nlayer),
+     $        b007(nlayer), b008(nlayer), b009(nlayer)
+      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
+     $        c004(nlayer), c005(nlayer), c006(nlayer), 
+     $        c007(nlayer), c008(nlayer), c009(nlayer), 
+     $        c010(nlayer), c011(nlayer), c012(nlayer),
+     $        c013(nlayer), c014(nlayer), c015(nlayer), 
+     $        c016(nlayer), c017(nlayer), c018(nlayer)
+      real :: d001(nlayer), d002(nlayer), d003(nlayer)
+      real :: e001(nlayer), e002(nlayer), e003(nlayer)
+      real :: h001(nlayer), h002(nlayer), h003(nlayer),
+     $        h004(nlayer), h005(nlayer)
+      real :: t001(nlayer), t002(nlayer), t003(nlayer)
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     stimestep  : standard timestep for the chemistry (s)       c
@@ -91,22 +93,23 @@
 c     phychemrat : ratio physical/chemical timestep              c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
+
       stimestep = 600. ! standard value : 10 mn
-c
+
       phychemrat = nint(ptimestep/stimestep)
-c
+
       ctimestep = ptimestep/real(phychemrat)
-c
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     initialisation of chemical species and families            c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call gcmtochim(zycol, lswitch, nesp, rm) 
-c
+
+      call gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm) 
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     compute reaction rates                                     c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call chemrates(lswitch, dens, press, temp,
+
+      call chemrates(nlayer,
+     $               lswitch, dens, press, temp,
      $               surfdust1d, surfice1d,
      $               a001, a002, a003,
@@ -120,32 +123,32 @@
      $               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, 
+
+      call phot(nlayer, 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, 
+         call chimie(nlayer, lswitch, nesp, rm, j, dens, ctimestep, 
      $               press, temp, sza, dist_sol,
      $               a001, a002, a003,
@@ -160,17 +163,18 @@
      $               t001, t002, t003)
       end do
-c
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     save chemical species for the gcm                          c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      call chimtogcm(zycol, lswitch, nesp, rm) 
-c
+
+      call chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm) 
+
       return
       end
-c
+
 c*****************************************************************
-c
-      subroutine chimie(lswitch, nesp, rm, j, dens, dt, 
+
+      subroutine chimie(nlayer,
+     $                  lswitch, nesp, rm, j, dens, dt, 
      $                  press, t, sza, dist_sol,
      $                  a001, a002, a003,
@@ -184,90 +188,88 @@
      $                  h001, h002, h003, h004, h005,
      $                  t001, t002, t003)
-c
+
 c*****************************************************************
-c
+
       implicit none 
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "chimiedata.h" 
 #include "callkeys.h" 
-c
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     inputs:
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer :: lswitch            ! interface level between chemistries
+      integer :: nesp               ! number of species
+
+      real :: dens(nlayer)          ! density (cm-3)
+      real :: dt                    ! chemistry timestep (s) 
+      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
+      real :: press(nlayer)         ! pressure (hpa)
+      real :: t(nlayer)             ! temperature (k)
+      real :: sza                   ! solar zenith angle (deg)
+      real :: dist_sol              ! sun distance (au) 
+
+c     reaction rates
+
+      real :: a001(nlayer), a002(nlayer), a003(nlayer)
+      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
+     $        b004(nlayer), b005(nlayer), b006(nlayer),
+     $        b007(nlayer), b008(nlayer), b009(nlayer)
+      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
+     $        c004(nlayer), c005(nlayer), c006(nlayer), 
+     $        c007(nlayer), c008(nlayer), c009(nlayer), 
+     $        c010(nlayer), c011(nlayer), c012(nlayer),
+     $        c013(nlayer), c014(nlayer), c015(nlayer), 
+     $        c016(nlayer), c017(nlayer), c018(nlayer)
+      real :: d001(nlayer), d002(nlayer), d003(nlayer)
+      real :: e001(nlayer), e002(nlayer), e003(nlayer)
+      real :: h001(nlayer), h002(nlayer), h003(nlayer),
+     $        h004(nlayer), h005(nlayer)
+      real :: t001(nlayer), t002(nlayer), t003(nlayer)
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     input/output:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-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
+
+      real :: rm(nlayer,nesp)  ! volume mixing ratios 
+
 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
+
+      real, parameter :: hetero_ice  = 1. ! switch for het. chem. on ice clouds
+      real, parameter :: hetero_dust = 0. ! switch for het. chem. on dust
+
+      integer :: iesp, iter, l
+      integer :: i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
+     $           i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
+      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
+     $           j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
+     $           j_ho2, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
+     $           j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
+     $           j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
+     $           j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
+     $           j_hno3, j_hno4 
+
+      integer, parameter :: niter = 5        ! iterations in the chemical scheme
+
+      real :: cc0(nlayer,nesp)        ! initial number density (cm-3)
+      real :: cc(nlayer,nesp)         ! final number density (cm-3)
+      real :: nox(nlayer)             ! nox number density (cm-3)
+      real :: no(nlayer)              ! no  number density (cm-3)
+      real :: no2(nlayer)             ! no2 number density (cm-3)
+      real :: production(nlayer,nesp) ! production rate
+      real :: loss(nlayer,nesp)       ! loss rate
+
+      real :: ro_o3, rh_ho2, roh_ho2, rno2_no
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     tracer numbering in the chemistry
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
+
       i_co2    =  1
       i_co     =  2
@@ -286,9 +288,9 @@
       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)
@@ -320,9 +322,9 @@
       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
@@ -331,5 +333,5 @@
          end do
       end do
-c
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     co2 and nox number densities (cm-3)    
@@ -337,15 +339,15 @@
 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
@@ -353,25 +355,25 @@
 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
@@ -381,7 +383,7 @@
 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)
@@ -400,9 +402,9 @@
      $            /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
@@ -428,26 +430,26 @@
      $            + 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
@@ -462,13 +464,13 @@
 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)
@@ -480,9 +482,9 @@
      $                   + 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)
@@ -492,21 +494,21 @@
      $              )
      $              /(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)
@@ -518,5 +520,5 @@
      $                   + 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)
@@ -529,38 +531,37 @@
      $                   + 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)
@@ -570,42 +571,42 @@
      $                   + 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)
@@ -625,24 +626,24 @@
      $                + 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)
@@ -652,13 +653,13 @@
      $                 + 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)
@@ -667,5 +668,5 @@
 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)
@@ -673,9 +674,9 @@
      $                  + 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)
@@ -685,5 +686,5 @@
      $                 + 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)
@@ -697,13 +698,13 @@
      $                 + 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)
@@ -714,13 +715,13 @@
      $                 + 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)
@@ -741,7 +742,7 @@
      $              /(1. + loss(l,i_ch4)*dt)
       end do
-c
+
 c     ox species
-c
+
       if (sza .le. 95.) then
          do l = 1,lswitch-1
@@ -758,15 +759,15 @@
          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
@@ -774,78 +775,80 @@
          end do
       end do
-c
+
       return
       end
-c
+
 c*****************************************************************
-c
-      subroutine phot(lswitch, press, temp, sza, tauref, dist_sol, 
+
+      subroutine phot(nlayer,
+     $                lswitch, press, temp, sza, tauref, dist_sol, 
      $                rmco2, rmo3, j)
-c
+
 c*****************************************************************
-c
+
       implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#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
+       
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer :: lswitch            ! interface level between chemistries
+      real :: press(nlayer)         ! pressure (hPa)
+      real :: temp(nlayer)          ! temperature (K)
+      real :: sza                   ! solar zenith angle (deg)
+      real :: tauref                ! optical depth at 7 hpa
+      real :: dist_sol              ! sun distance (AU) 
+      real :: rmco2(nlayer)         ! co2 mixing ratio
+      real :: rmo3(nlayer)          ! ozone mixing ratio
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     output:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real j(nlayermx,nd)              ! interpolated photolysis rates (s-1) 
-c
+
+      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
+
 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
+
+      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp, 
+     $           iozo, isza, itau, it, l
+
+      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
+      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
+      real :: poids(2,2,2,2,2)            ! 5D interpolation weights 
+      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
+      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
+      real :: cinf, csup, cicol, 
+     $        ciozo, cisza, citemp, citau
+      real :: colo3min, dp, coef
+      real :: ratio_o3(nlayer)
+      real :: tau
+
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     day/night criterion
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-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
@@ -857,12 +860,12 @@
       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
@@ -874,20 +877,20 @@
       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.
@@ -900,5 +903,5 @@
      $              *6.022e22*dp/(mugaz*g)
       end do
-c
+
 c     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
 c
@@ -906,5 +909,5 @@
 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
@@ -913,11 +916,11 @@
          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
@@ -931,10 +934,10 @@
       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
@@ -946,13 +949,13 @@
       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
@@ -964,7 +967,7 @@
          end do
  30      continue
-c
+
 cc    interpolation in ozone column
-c
+
          indozo = nozo - 1
          do iozo = 1,nozo
@@ -976,9 +979,9 @@
          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)
@@ -1013,5 +1016,5 @@
          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
@@ -1046,7 +1049,7 @@
          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) =
@@ -1083,5 +1086,5 @@
      $   + 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)
@@ -1117,23 +1120,23 @@
      $     *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
@@ -1141,47 +1144,50 @@
          end do
       end do
-c
+
       end if
-c
+
       return
       end
-c
+
 c*****************************************************************
-c
-      subroutine gcmtochim(zycol, lswitch, nesp, rm)
-c
+
+      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
+
 c*****************************************************************
-c
-      use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d,
+
+      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
      &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
      &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
      &                      igcm_ch4
+
       implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#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
+      
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer, intent(in) :: nq     ! number of tracers
+      real :: zycol(nlayer,nq)      ! species volume mixing ratio in the gcm
+
+      integer :: nesp               ! number of species in the chemistry
+      integer :: lswitch            ! interface level between chemistries
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     outputs:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      real rm(nlayermx,nesp)   ! species volume mixing ratio 
-c      
+
+      real :: rm(nlayer,nesp)   ! species volume mixing ratio 
+      
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     local:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer l,iq
+
+      integer :: l, iq
       
 c     tracer indexes in the chemistry:
@@ -1203,9 +1209,9 @@
       integer,parameter :: i_hox  = 15
       integer,parameter :: i_ox   = 16
-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)
@@ -1233,9 +1239,9 @@
          end do
       end if
-c
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     initialise chemical families                     c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
+
       do l = 1,lswitch-1
          rm(l,i_hox) = rm(l,i_h) 
@@ -1245,45 +1251,48 @@
      $               + rm(l,i_o3)
       end do
-c
+
       return
       end
-c
+
 c*****************************************************************
 c
-      subroutine chimtogcm(zycol, lswitch, nesp, rm)
+      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
 c
 c*****************************************************************
 c
-      use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d,
+      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
      &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
      &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
      &                      igcm_ch4
+
       implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#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
+
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer, intent(in) :: nq     ! number of tracers
+      integer :: nesp               ! number of species in the chemistry
+      integer :: lswitch            ! interface level between chemistries
+      
+      real :: rm(nlayer,nesp)       ! species volume mixing ratio 
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     output:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c      
-      real zycol(nlayermx,nqmx)  ! species volume mixing ratio in the gcm
-c
+
+      real :: zycol(nlayer,nq)  ! species volume mixing ratio in the gcm
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     local:
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer l,iq
+
+      integer l, iq
       
 c     tracer indexes in the chemistry:
@@ -1305,9 +1314,9 @@
       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) 
@@ -1331,11 +1340,12 @@
          end do
       end if
-c
+
       return
       end
-c
+
 c*****************************************************************
-c
-      subroutine chemrates(lswitch, dens, press, t,
+
+      subroutine chemrates(nlayer,
+     $                     lswitch, dens, press, t,
      $                     surfdust1d, surfice1d,
      $                     a001, a002, a003,
@@ -1349,66 +1359,66 @@
      $                     h001, h002, h003, h004, h005,
      $                     t001, t002, t003, tau)
-c
+
 c*****************************************************************
-c
+
       implicit none
-c
-#include "dimensions.h"
-#include "dimphys.h"
+
+!#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
+
+      integer, intent(in) :: nlayer   ! number of atmospheric layers
+      integer :: lswitch              ! interface level between chemistries
+
+      real :: dens(nlayer)            ! density (cm-3)
+      real :: press(nlayer)           ! pressure (hpa)
+      real :: t(nlayer)               ! temperature (k)
+      real :: surfdust1d(nlayer)      ! dust surface area (cm^2/cm^3)
+      real :: surfice1d(nlayer)       ! ice surface area (cm^2/cm^3)
+      real :: tau                     ! dust opacity at 7 hpa
+
+      real, parameter :: tribo   = 0. ! switch for triboelectricity
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     outputs:                                                   c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-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
+
+      real :: a001(nlayer), a002(nlayer), a003(nlayer)
+      real :: b001(nlayer), b002(nlayer), b003(nlayer), 
+     $        b004(nlayer), b005(nlayer), b006(nlayer),
+     $        b007(nlayer), b008(nlayer), b009(nlayer)
+      real :: c001(nlayer), c002(nlayer), c003(nlayer), 
+     $        c004(nlayer), c005(nlayer), c006(nlayer), 
+     $        c007(nlayer), c008(nlayer), c009(nlayer), 
+     $        c010(nlayer), c011(nlayer), c012(nlayer),
+     $        c013(nlayer), c014(nlayer), c015(nlayer), 
+     $        c016(nlayer), c017(nlayer), c018(nlayer)
+      real :: d001(nlayer), d002(nlayer), d003(nlayer)
+      real :: e001(nlayer), e002(nlayer), e003(nlayer)
+      real :: h001(nlayer), h002(nlayer), h003(nlayer),
+     $        h004(nlayer), h005(nlayer)
+      real :: t001(nlayer), t002(nlayer), t003(nlayer)
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     local:                                                     c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-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
+
+      real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
+      real :: ef, efmax, lossh2o, lossch4, lossco2
+
+      integer :: l
+      real :: k1a, k1b, k1a0, k1b0, k1ainf
+      real :: x, y, fc, fx
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     compute reaction rates
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
+
       do l = 1,lswitch-1
-c
+
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c        oxygen compounds
