SUBROUTINE chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
     $                       zlocal,zenit,ptimestep,zday)

       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing chemical variations in the thermosphere
c
c   author:  MAC July 2003
c   ------
c   modifications:
c   -------------
c     Ehouarn Sept 2008: added handling of tracers by their names
c=======================================================================
c
c    0.  Declarations :
c    ------------------
c
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "comdiurn.h"
#include "param.h"
#include "param_v3.h"
#include "chimiedata.h"
#include "conc.h"
#include"tracer.h"
c-----------------------------------------------------------------------
c    Input/Output
c    ------------
      INTEGER lswitch,ig

      REAL zday,zycol(nlayermx,nqmx)
      REAL ptimestep
      real zenit
      real ztemp(nlayermx)
      real zdens(nlayermx)
      real zpress(nlayermx)			! in mbar
c
c    Local variables :
c    -----------------
      INTEGER nlayer,l,nesptherm,iq
      parameter (nesptherm = 11)
      real tmean
      real aux1(nlayermx),aux2(nlayermx)
      real zlocal(nlayermx)
      real rm(nlayermx,nesptherm) 		!number density (cm-3)

      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 thermospheric chemistry:
      integer,parameter :: i_co2=1
      integer,parameter :: i_o2=2
      integer,parameter :: i_o=3
      integer,parameter :: i_co=4
      integer,parameter :: i_h=5
      integer,parameter :: i_oh=6
      integer,parameter :: i_ho2=7
      integer,parameter :: i_h2=8
      integer,parameter :: i_h2o=9
      integer,parameter :: i_h2o2=10
      integer,parameter :: i_o1d=11


! Initializations at first call
      if (firstcall) then
c        call param_read
        ! get the indexes of the tracers we'll need
        g_co2=igcm_co2
        if (g_co2.eq.0) then
          write(*,*) "chemthermos: Error; no CO2 tracer !!!"
          stop
        endif
        g_co=igcm_co
        if (g_co.eq.0) then
          write(*,*) "chemthermos: Error; no CO tracer !!!"
          stop
        endif
        g_o=igcm_o
        if (g_o.eq.0) then
          write(*,*) "chemthermos: Error; no O tracer !!!"
          stop
        endif
        g_o1d=igcm_o1d
        if (g_o1d.eq.0) then
          write(*,*) "chemthermos: Error; no O1D tracer !!!"
          stop
        endif
        g_o2=igcm_o2
        if (g_o2.eq.0) then
          write(*,*) "chemthermos: Error; no O2 tracer !!!"
          stop
        endif
        g_h=igcm_h
        if (g_h.eq.0) then
          write(*,*) "chemthermos: Error; no H tracer !!!"
          stop
        endif
        g_h2=igcm_h2
        if (g_h2.eq.0) then
          write(*,*) "chemthermos: Error; no H2 tracer !!!"
          stop
        endif
        g_oh=igcm_oh
        if (g_oh.eq.0) then
          write(*,*) "chemthermos: Error; no OH tracer !!!"
          stop
        endif
        g_ho2=igcm_ho2
        if (g_ho2.eq.0) then
          write(*,*) "chemthermos: Error; no HO2 tracer !!!"
          stop
        endif
        g_h2o2=igcm_h2o2
        if (g_h2o2.eq.0) then
          write(*,*) "chemthermos: Error; no H2O2 tracer !!!"
          stop
        endif
        g_h2o_vap=igcm_h2o_vap
        if (g_h2o_vap.eq.0) then
          write(*,*) "chemthermos: Error; no water vapor tracer !!!"
          stop
        endif

        firstcall= .false.
      endif ! of if (firstcall)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      nlayer=nlayermx
c      zlocal(1)=0.00625
c      do l=2,nlayer
c        tmean=ztemp(l)
c        if(ztemp(l).ne.ztemp(l-1))
c     &    tmean=(ztemp(l)-ztemp(l-1))/log(ztemp(l)/ztemp(l-1))
c        zlocal(l)= zlocal(l-1)-log(zpress(l)/zpress(l-1))
c     &                         *Rnew(ig,l-1)*tmean/g/1000.
c      enddo

      do l=1,nlayer
        aux1(l)=0.
        aux2(l)=0.
        rm(l,i_co2)  = zycol(l,g_co2)    *zdens(l)  
        rm(l,i_co)   = zycol(l,g_co)     *zdens(l) 
        rm(l,i_o)    = zycol(l,g_o)      *zdens(l)
        rm(l,i_o1d)  = zycol(l,g_o1d)    *zdens(l) 
        rm(l,i_o2)   = zycol(l,g_o2)     *zdens(l)
        rm(l,i_h)    = zycol(l,g_h)      *zdens(l) 
        rm(l,i_h2)   = zycol(l,g_h2)     *zdens(l) 
        rm(l,i_oh)   = zycol(l,g_oh)     *zdens(l) 
        rm(l,i_ho2)  = zycol(l,g_ho2)    *zdens(l)
        rm(l,i_h2o2) = zycol(l,g_h2o2)   *zdens(l) 
        rm(l,i_h2o)  = zycol(l,g_h2o_vap)*zdens(l) 
      enddo
         
      call flujo(solarcondate+zday/365.)

      call jthermcalc
     $ (rm(1,i_co2),rm(1,i_o2),rm(1,i_o),rm(1,i_h2),rm(1,i_h2o),
     &  rm(1,i_h2o2),aux1,aux2,ztemp,nlayermx,zlocal,
     &  solarcondate+zday/365.,zenit)

      call paramfoto(lswitch,zdens,ztemp,ptimestep/3600.,zenit,
     &             nlayer,rm(1,i_co2),rm(1,i_o2),rm(1,i_o),
     &             rm(1,i_co),rm(1,i_h),rm(1,i_oh),rm(1,i_ho2),
     &             rm(1,i_h2),rm(1,i_h2o),rm(1,i_h2o2),rm(1,i_o1d))

      do l=lswitch,nlayer
        zycol(l,g_co2)     = max(rm(l,i_co2)  / zdens(l) , 1.e-30)
        zycol(l,g_co)      = max(rm(l,i_co)   / zdens(l) , 1.e-30)
        zycol(l,g_o2)      = max(rm(l,i_o2)   / zdens(l) , 1.e-30)
        zycol(l,g_h2)      = max(rm(l,i_h2)   / zdens(l) , 1.e-30)
        zycol(l,g_h)       = max(rm(l,i_h)    / zdens(l) , 1.e-30)
        zycol(l,g_oh)      = max(rm(l,i_oh)   / zdens(l) , 1.e-30)
        zycol(l,g_ho2)     = max(rm(l,i_ho2)  / zdens(l) , 1.e-30)
        zycol(l,g_h2o_vap) = max(rm(l,i_h2o)  / zdens(l) , 1.e-30)
        zycol(l,g_h2o2)    = max(rm(l,i_h2o2) / zdens(l) , 1.e-30)
        zycol(l,g_o1d)     = max(rm(l,i_o1d)  / zdens(l) , 1.e-30)
        zycol(l,g_o)       = max(rm(l,i_o)    / zdens(l) , 1.e-30)
       enddo	!nlayer


      return
      end