subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
     $                   zzlay,zday,pq,pdq,rice,
     &                   dqchim,dqschim,dqcloud,dqscloud)
c
      implicit none
c
c=======================================================================
c
c   subject:
c   --------
c
c  Prepare the call for the photochemical module, and send back the
c  tendencies from photochemistry in the chemical species mass mixing ratios
c
c   Author:   Sebastien Lebonnois (08/11/2002)
c   -------
c    update 12/06/2003 for water ice clouds and compatibility with dust
c    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
c    update 03/05/2005 cosmetic changes (Franck Lefevre)
c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
c
c   Arguments:
c   ----------
c
c  Input:
c
c    ptimestep                  timestep (s)
c    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
c    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
c    pt(ngridmx,nlayermx)       Temperature (K)
c    pdt(ngridmx,nlayermx)      Temperature 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
c  Output:
c
c    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
c    dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
c
c=======================================================================

c    Declarations :
c    --------------

#include "dimensions.h"
#include "dimphys.h"
#include "chimiedata.h"
#include "tracer.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "conc.h"

c Arguments :
c -----------

c   inputs:
c   -------

      real    ptimestep
      real    pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
      real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
      real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
      real    pt(ngridmx,nlayermx)       ! temperature
      real    pdt(ngridmx,nlayermx)      ! temperature tendency
      real    dist_sol                   ! distance of the sun (AU)
      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)

c   outputs:
c   --------

      real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
      real dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
      real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
      real dqscloud(ngridmx,nqmx)        ! tendencies on qsurf 

c Local variables :
c -----------------

      character*20 str20
      integer  ig,l,i,iq
      integer  foundswitch, lswitch
      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
c
c  for each column of atmosphere:
c
      real zpress(nlayermx)       !  Pressure (mbar)
      real zdens(nlayermx)        !  Density  (cm-3)
      real ztemp(nlayermx)        !  Temperature (K)
      real zlocal(nlayermx)       !  Altitude (km)
      real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
      real szacol                 !  Solar zenith angle
      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
      parameter (output=.true.)   ! see at end of routine

      logical,save :: firstcall=.true.
      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
      integer,save :: i_o=0
      integer,save :: i_o1d=0
      integer,save :: i_o2=0
      integer,save :: i_o3=0
      integer,save :: i_h=0
      integer,save :: i_h2=0
      integer,save :: i_oh=0
      integer,save :: i_ho2=0
      integer,save :: i_h2o2=0
      integer,save :: i_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=======================================================================
c     initialization of the chemistry (first call only)
c=======================================================================
c
      if (firstcall) then
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
         end if

         ! find index of chemical tracers to use
         nbq=0 ! to count number of tracers
         i_co2=igcm_co2
         if (i_co2.eq.0) then
           write(*,*) "calchim: Error; no CO2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_co2
         endif
         i_co=igcm_co
         if (i_co.eq.0) then
           write(*,*) "calchim: Error; no CO tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_co
         endif
         i_o=igcm_o
         if (i_o.eq.0) then
           write(*,*) "calchim: Error; no O tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o
         endif
         i_o1d=igcm_o1d
         if (i_o1d.eq.0) then
           write(*,*) "calchim: Error; no O1D tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o1d
         endif
         i_o2=igcm_o2
         if (i_o2.eq.0) then
           write(*,*) "calchim: Error; no O2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o2
         endif
         i_o3=igcm_o3
         if (i_o3.eq.0) then
           write(*,*) "calchim: Error; no O3 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o3
         endif
         i_h=igcm_h
         if (i_h.eq.0) then
           write(*,*) "calchim: Error; no H tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h
         endif
         i_h2=igcm_h2
         if (i_h2.eq.0) then
           write(*,*) "calchim: Error; no H2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2
         endif
         i_oh=igcm_oh
         if (i_oh.eq.0) then
           write(*,*) "calchim: Error; no OH tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_oh
         endif
         i_ho2=igcm_ho2
         if (i_ho2.eq.0) then
           write(*,*) "calchim: Error; no HO2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ho2
         endif
         i_h2o2=igcm_h2o2
         if (i_h2o2.eq.0) then
           write(*,*) "calchim: Error; no H2O2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2o2
         endif
         i_n2=igcm_n2
         if (i_n2.eq.0) then
           write(*,*) "calchim: Error; no N2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_n2
         endif
         i_ar=igcm_ar
         if (i_ar.eq.0) then
           write(*,*) "calchim: Error; no AR tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ar
         endif
         i_ice=igcm_h2o_ice
         if (i_ice.eq.0) then
           write(*,*) "calchim: Error; no water ice tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ice
         endif
         i_h2o=igcm_h2o_vap
         if (i_h2o.eq.0) then
           write(*,*) "calchim: Error; no water vapor tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2o
         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
         
         firstcall = .false.
      end if ! if (firstcall)

! Initialize output tendencies to zero (to handle case of tracers which
! are not used in the chemistry (e.g. dust))
      dqchim(:,:,:)=0
      dqschim(:,:)=0


c
c=======================================================================
c     loop over grid
c=======================================================================
c
      do ig = 1,ngridmx
c
c     local updates
c 
         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
              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
            zpress(l) = pplay(ig,l)/100.
            ztemp(l)  = zt(ig,l)
            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
            zlocal(l) = zzlay(ig,l)/1000.
c
c     search for switch index between regions
c
            if (photochem .and. thermochem) then
               if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
                  lswitch = l
                  foundswitch=1
               end if
            end if
            if ( .not. photochem) then
               lswitch = 22
            end if
            if (.not.  thermochem) then
               lswitch = min(33,nlayermx+1)
            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
c
c=======================================================================
c     call chemical subroutine
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
c=======================================================================
c     tendencies
c=======================================================================
c
c     must be 0. for water ice:
c
         if (water) then
            do l = 1,nlayermx
!               dqchim(ig,l,nqmx-1) = 0.
               dqchim(ig,l,i_ice) = 0.
            end do
         end if
c
c     tendency for CO2 = - sum of others for lower atmosphere
c     tendency for O   = - sum of others for upper atmosphere
c
         do l = 1,nlayermx
            if (l .lt. lswitch) then
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                  if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then
                     dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)
     $                                - zq(ig,l,iq))/ptimestep
                  else if (iq.eq.i_co2) then
                     dqchim(ig,l,iq) = 0.
                  end if
                  dqschim(ig,iq) = 0.
               end do ! of do i=1,nbq
               
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                  if (iq.ne.i_co2) then
                     dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2) 
     $                                  - dqchim(ig,l,iq)
                  end if
               end do
             else if (l .ge. lswitch) then
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                   if ((iq.ne.i_o).and.(iq.ne.i_ice)) then
                      dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)
     $                                  /mmean(ig,l)
     $                                 - zq(ig,l,iq))/ptimestep
                   else if (iq.eq.i_o) then
!                      i_o = iq
                      dqchim(ig,l,iq) = 0.
                   end if
               enddo
               
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                   if (iq.ne.i_o) then
                      dqchim(ig,l,i_o) = dqchim(ig,l,i_o) 
     $                                 - dqchim(ig,l,iq)
                   end if
               end do
             end if ! of if (l.lt.lswitch) else if (l.ge.lswitch)
          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
          call perosat(ig,ptimestep,pplev,pplay,
     $                 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)
             jo3_3d(ig,l) = jo3(l)
          end do
c
c=======================================================================
c     end of loop over grid
c=======================================================================
c
      end do ! of do ig=1,ngridmx
c
c=======================================================================
c     write outputs
c=======================================================================
c
! value of parameter 'output' to trigger writting of outputs
! is set above at the declaration of the variable.

      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)
         end if ! of if (ngridmx.gt.1)
c
      endif ! of if (output)
c
      end