SUBROUTINE concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
                                             
      implicit none

c=======================================================================
c CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
c
c mmean(ngridmx,nlayermx)	amu
c cpnew(ngridmx,nlayermx)	J/kg/K
c rnew(ngridmx,nlayermx)	J/kg/K
c akknew(ngridmx,nlayermx)	coefficient of thermal concduction
c
c version: March 2011 - Franck Lefevre
c=======================================================================

c    Declarations
c    ------------
 
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "comdiurn.h"
#include "chimiedata.h"
#include "tracer.h"
#include "conc.h"

c    Input/Output
c    ------------

      real pplay(ngridmx,nlayermx)
      real pt(ngridmx,nlayermx)
      real pdt(ngridmx,nlayermx)
      real pq(ngridmx,nlayermx,nqmx)
      real pdq(ngridmx,nlayermx,nqmx)
      real ptimestep

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

      integer       :: l, ig, n, k
      integer, save :: gind(ncomp)
      real          :: ni(nqmx), ntot
      real          :: zq(ngridmx,nlayermx,ncomp)
      real          :: zt(ngridmx,nlayermx)
      real, save    :: aki(ncomp)
      real, save    :: cpi(ncomp)

      logical, save :: firstcall = .true.

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     tracer numbering for the thermal conduction and
c     specific heat coefficients
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      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_n2   = 13
      integer,parameter :: i_ar   = 14
      integer,parameter :: i_h2o  = 15

      if (firstcall) then

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c        initializations at first call:
c        fill local array of tracer indexes
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

         gind(i_co2)   =  igcm_co2     ! co2
         gind(i_co)    =  igcm_co      ! co
         gind(i_o)     =  igcm_o       ! o
         gind(i_o1d)   =  igcm_o1d     ! o1d
         gind(i_o2)    =  igcm_o2      ! o2
         gind(i_o3)    =  igcm_o3      ! o3
         gind(i_h)     =  igcm_h       ! h
         gind(i_h2)    =  igcm_h2      ! h2
         gind(i_oh)    =  igcm_oh      ! oh
         gind(i_ho2)   =  igcm_ho2     ! ho2
         gind(i_h2o2)  =  igcm_h2o2    ! h2o2
         gind(i_ch4)   =  igcm_ch4     ! ch4
         gind(i_n2)    =  igcm_n2      ! n2
         gind(i_ar)    =  igcm_ar      ! ar
         gind(i_h2o)   =  igcm_h2o_vap ! h2o

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c    Thermal conductivity and specific heat coefficients
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

         aki(i_co2)   = 3.072e-4
         aki(i_co)    = 4.87e-4
         aki(i_o)     = 7.59e-4
         aki(i_o1d)   = 7.59e-4    !?
         aki(i_o2)    = 5.68e-4
         aki(i_o3)    = 3.00e-4    !?
         aki(i_h)     = 0.0
         aki(i_h2)    = 36.314e-4
         aki(i_oh)    = 7.00e-4    !?
         aki(i_ho2)   = 0.0
         aki(i_h2o2)  = 0.0
         aki(i_ch4)   = 0.0        !?
         aki(i_n2)    = 5.6e-4
         aki(i_ar)    = 0.0        !?
         aki(i_h2o)   = 0.0

         cpi(i_co2)   = 0.834e3
         cpi(i_co)    = 1.034e3
         cpi(i_o)     = 1.3e3
         cpi(i_o1d)   = 1.3e3    !?
         cpi(i_o2)    = 0.9194e3
         cpi(i_o3)    = 0.800e3  !?
         cpi(i_h)     = 20.780e3
         cpi(i_h2)    = 14.266e3
         cpi(i_oh)    = 1.045e3
         cpi(i_ho2)   = 1.065e3  !?
         cpi(i_h2o2)  = 1.000e3  !?
         cpi(i_ch4)   = 1.000e3  !?
         cpi(i_n2)    = 1.034e3
         cpi(i_ar)    = 1.000e3  !?
         cpi(i_h2o)   = 1.870e3

         firstcall=.false.

      end if ! of if (firstcall)
c
c     initializations
c
      mmean(:,:)  = 0.
      cpnew(:,:)  = 0.
      akknew(:,:) = 0.
c
c     update temperature
c
      do l = 1,nlayermx
         do ig = 1,ngridmx
            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
         end do
      end do
c
c     update tracers
c
      do l = 1,nlayermx
         do ig = 1,ngridmx
            do n = 1,ncomp
               zq(ig,l,n) = max(1.e-30, pq(ig,l,gind(n))
     $                                + pdq(ig,l,gind(n))*ptimestep)
            end do
         end do
      end do
c
c     mmean : mean molecular mass
c     rnew  : specific gas constant
c
      do l = 1,nlayermx
         do ig = 1,ngridmx
            do n = 1, ncomp
               mmean(ig,l) = mmean(ig,l) + zq(ig,l,n)/mmol(gind(n))
            end do
            mmean(ig,l) = 1./mmean(ig,l)
            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg/K		
         end do
      end do
c
c     cpnew  : specicic heat
c     akknew : thermal conductivity cofficient
c      
      do l = 1,nlayermx
         do ig = 1,ngridmx
            ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
            do n = 1,ncomp
               ni(n) = ntot*zq(ig,l,n)*mmean(ig,l)/mmol(gind(n))
               cpnew(ig,l) = cpnew(ig,l) + ni(n)*cpi(n)
               akknew(ig,l) = akknew(ig,l) + ni(n)*aki(n)
            end do 
            cpnew(ig,l) = cpnew(ig,l)/ntot
            akknew(ig,l) = akknew(ig,l)/ntot
         end do
c        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
      end do

      return
      end