source: trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F @ 1314

      SUBROUTINE concentrations2(pplay,t_seri,pdt,co2vmr_gcm, n2vmr_gcm,
     $   covmr_gcm,o3pvmr_gcm,nvmr_gcm,ptimestep)

      use dimphy
      use conc,  only: mmean, Akknew, rnew, cpnew
      use cpdet_mod, only: cpdet                       
      implicit none

!=======================================================================
! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
!
! mmean(klon,klev)      amu
! cpnew(klon,klev)      J/kg/K
! rnew(klon,klev)       J/kg/K
! akknew(klon,klev)     coefficient of thermal conduction
!
! version: April 2012 - Franck Lefevre
!=======================================================================

!     declarations
 
#include "dimensions.h"
#include "YOMCST.h"
#include "clesphys.h"
c#include "comdiurn.h"
c#include "chimiedata.h"
c#include "tracer.h"
#include "mmol.h"

!     input/output

      real pplay(klon,klev)
c      real pt(klon,klev)
      real t_seri(klon, klev)
      real pdt(klon,klev)
      real co2vmr_gcm(klon,klev), n2vmr_gcm(klon,klev)
      real covmr_gcm(klon,klev)
      real o3pvmr_gcm(klon,klev),nvmr_gcm(klon,klev)
c      real pq(klon,klev,nqmx)
c      real pdq(klon,klev,nqmx)
      real ptimestep

!     local variables

      integer       :: i, l, ig, iq
      real          :: ntot
      real          :: zt(klon, klev)
      real, save    :: akico2,akio,akin2,akico
      real, save    :: akin
      real, save    :: cpico2,cpico,cpio,cpin2
      real, save    :: cpio2

      logical, save :: firstcall = .true.

      if (firstcall) then

!        initialize thermal conductivity and specific heat coefficients
!        values are taken from the literature [J/kg K]

 !          co2 
            akico2 = 3.072e-4
            cpico2 = 0.834e3

 !       co
            akico = 4.87e-4
            cpico = 1.034e3

 !        o
            akio = 7.59e-4
            cpio = 1.3e3

 !       n 
 !           akin = 0.0
 !           cpin = 0.0

  !     n2 
            akin2 = 5.6e-4
            cpin2 = 1.034e3

         firstcall = .false.
      end if ! if (firstcall)

!     update temperature

      do l = 1,klev
         do ig = 1,klon
c            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
            zt(ig,l) = t_seri(ig,l) 

         end do
      end do

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

      mmean(:,:)  = 0.

      do l = 1,klev
         do ig = 1,klon
c            do i = 1,nbq
c               iq = niq(i) 
c               mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
c            end do
            mmean(ig,l) = RMD
            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg K            
         end do
      end do

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

      do l = 1,klev
          do ig = 1,klon

            ntot = pplay(ig,l)/(RKBOL*zt(ig,l))*1.e-6  ! in #/cm3 

            cpnew(ig,l) = ntot*o3pvmr_gcm(ig,l)*cpio 
     $              +ntot* co2vmr_gcm(ig,l)*cpdet(zt(ig,l))             
     $         + ntot*n2vmr_gcm(ig,l)*cpin2 + ntot*covmr_gcm(ig,l)*cpico    


            akknew(ig,l) = ntot*o3pvmr_gcm(ig,l)*akio + 
     $               ntot*co2vmr_gcm(ig,l)*akico2 +
     $           ntot*n2vmr_gcm(ig,l)*akin2 + ntot*covmr_gcm(ig,l)*akico          


            cpnew(ig,l) = cpnew(ig,l)/ntot
            akknew(ig,l)= akknew(ig,l)/ntot

c        print*, '--- concentrations ---'
c        print*, l, cpnew(1,l),  rnew(1,l), akknew(1,l)
          end do 
       end do
c       STOP

      return
      end