SUBROUTINE concentrations(ngrid,nlayer,nq,
     &                          pplay,pt,pdt,pq,pdq,ptimestep)
                                             
      use tracer_h, only: mmol, noms, aki, cpi, nesp

      use conc_mod, only: mmean, akknew, rnew, cpnew
      USE comcstfi_mod                   
      use callkeys_mod
      use chimiedata_h
      implicit none

!=======================================================================
! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
!
! mmean(ngrid,nlayer)	amu
! cpnew(ngrid,nlayer)	J/kg/K
! rnew(ngrid,nlayer)	J/kg/K
! akknew(ngrid,nlayer)	coefficient of thermal concduction
!
! version: April 2012 - Franck Lefevre
! update 06/03/2021 cpi/aki input (Yassin Jaziri)
!=======================================================================

!     input/output

      integer,intent(in) :: ngrid  ! number of atmospheric columns
      integer,intent(in) :: nlayer ! number of atmospheric layers
      integer,intent(in) :: nq     ! number of tracers
      real,   intent(in) :: pplay(ngrid,nlayer)
      real,   intent(in) :: pt(ngrid,nlayer)
      real,   intent(in) :: pdt(ngrid,nlayer)
      real,   intent(in) :: pq(ngrid,nlayer,nq)
      real,   intent(in) :: pdq(ngrid,nlayer,nq)
      real,   intent(in) :: ptimestep

!     local variables

      integer       :: l, ig, iq
      real          :: ni(nq), ntot
      real          :: zq(ngrid, nlayer, nq)
      real          :: zt(ngrid, nlayer)

!     update temperature

      do l = 1,nlayer
         do ig = 1,ngrid
            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
         end do
      end do

!     update tracers

      do l = 1,nlayer
         do ig = 1,ngrid
            do iq = 1,nq
               zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)
     $                                 + pdq(ig,l,iq)*ptimestep)
            end do
         end do
      end do

!     mmean : mean molecular mass
!     rnew  : specific gas constant

      mmean(:,:)  = 0.
      do l = 1,nlayer
         do ig = 1,ngrid
            do iq = 1,nq
               if (mmol(iq).ne.0.) then
                  mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
               end if
            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


!     cpnew  : specicic heat
!     akknew : thermal conductivity cofficient
      cpnew(:,:)  = 0.
      akknew(:,:) = 0.

      do l = 1,nlayer
         do ig = 1,ngrid
            ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
            do iq = 1,nq
               ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
               cpnew(ig,l)  = cpnew(ig,l) + ni(iq)*cpi(iq)
               akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(iq)
            end do 
            cpnew(ig,l)  = cpnew(ig,l)/ntot
            akknew(ig,l) = akknew(ig,l)/ntot
         end do
      end do

      return
      end