Index: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 358)
+++ trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 370)
@@ -1,18 +1,20 @@
       SUBROUTINE concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
                                              
-       IMPLICIT NONE
+      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 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=======================================================================
-c    0.  Declarations :
-c    ------------------
-c
+c    Declarations
+c    ------------
+ 
 #include "dimensions.h"
 #include "dimphys.h"
@@ -24,34 +26,32 @@
 #include "conc.h"
 
-c-----------------------------------------------------------------------
 c    Input/Output
 c    ------------
 
-      REAL pplay(ngridmx,nlayermx)
-      REAL pt(ngridmx,nlayermx)
-      REAL pdt(ngridmx,nlayermx)
+      real pplay(ngridmx,nlayermx)
+      real pt(ngridmx,nlayermx)
+      real pdt(ngridmx,nlayermx)
       real pq(ngridmx,nlayermx,nqmx)
-      REAL pdq(ngridmx,nlayermx,nqmx)
-      REAL ptimestep
+      real pdq(ngridmx,nlayermx,nqmx)
+      real ptimestep
 
-c    Local variables :
-c    -----------------
-      INTEGER,SAVE :: ngrid,nlayer,nq
-      INTEGER iq,l,ig,ll,n,k
-      integer,save :: gind(ncomptot)
-      real ni(ncomptot)
-      real nt, ntot
-      real q2(ngridmx,nlayermx,ncomptot)
-      real zt(ngridmx,nlayermx)
-      real q2tot(ngridmx,nlayermx)
-      real,save :: aki(ncomptot)
-      real,save :: cpi(ncomptot)
+c    Local variables
+c    ---------------
 
-      logical,save :: firstcall=.true.
+      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
-c
+
       integer,parameter :: i_co2  = 1
       integer,parameter :: i_co   = 2
@@ -65,202 +65,127 @@
       integer,parameter :: i_ho2  = 10
       integer,parameter :: i_h2o2 = 11
-      integer,parameter :: i_n2   = 12
-      integer,parameter :: i_ar   = 13
-      integer,parameter :: i_h2o  = 14
+      integer,parameter :: i_ch4  = 12
+      integer,parameter :: i_n2   = 13
+      integer,parameter :: i_ar   = 14
+      integer,parameter :: i_h2o  = 15
 
-! 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_ar=0
-      integer,save :: g_h2o_vap=0
-
-! Initializations at first call
       if (firstcall) then
-        ! identify the indexes of the tracers we'll need
-        g_co2=igcm_co2
-        if (g_co2.eq.0) then
-          write(*,*) "concentrations: Error; no CO2 tracer !!!"
-          stop
-        endif
-        g_co=igcm_co
-        if (g_co.eq.0) then
-          write(*,*) "concentrations: Error; no CO tracer !!!"
-          stop
-        endif
-        g_o=igcm_o
-        if (g_o.eq.0) then
-          write(*,*) "concentrations: Error; no O tracer !!!"
-          stop
-        endif
-        g_o1d=igcm_o1d
-        if (g_o1d.eq.0) then
-          write(*,*) "concentrations: Error; no O1D tracer !!!"
-          stop
-        endif
-        g_o2=igcm_o2
-        if (g_o2.eq.0) then
-          write(*,*) "concentrations: Error; no O2 tracer !!!"
-          stop
-        endif
-        g_o3=igcm_o3
-        if (g_o3.eq.0) then
-          write(*,*) "concentrations: Error; no O3 tracer !!!"
-          stop
-        endif
-        g_h=igcm_h
-        if (g_h.eq.0) then
-          write(*,*) "concentrations: Error; no H tracer !!!"
-          stop
-        endif
-        g_h2=igcm_h2
-        if (g_h2.eq.0) then
-          write(*,*) "concentrations: Error; no H2 tracer !!!"
-          stop
-        endif
-        g_oh=igcm_oh
-        if (g_oh.eq.0) then
-          write(*,*) "concentrations: Error; no OH tracer !!!"
-          stop
-        endif
-        g_ho2=igcm_ho2
-        if (g_ho2.eq.0) then
-          write(*,*) "concentrations: Error; no HO2 tracer !!!"
-          stop
-        endif
-        g_h2o2=igcm_h2o2
-        if (g_h2o2.eq.0) then
-          write(*,*) "concentrations: Error; no H2O2 tracer !!!"
-          stop
-        endif
-        g_n2=igcm_n2
-        if (g_n2.eq.0) then
-          write(*,*) "concentrations: Error; no N2 tracer !!!"
-          stop
-        endif
-        g_ar=igcm_ar
-        if (g_ar.eq.0) then
-          write(*,*) "concentrations: Error; no AR tracer !!!"
-          stop
-        endif
-        g_h2o_vap=igcm_h2o_vap
-        if (g_h2o_vap.eq.0) then
-          write(*,*) "concentrations: Error; no water vapor tracer !!!"
-          stop
-        endif
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     fill local array of tracer indexes
+c        initializations at first call:
+c        fill local array of tracer indexes
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-        gind(i_co2)   =  g_co2			! co2
-        gind(i_co)    =  g_co			! co
-        gind(i_o)     =  g_o			! o
-        gind(i_o1d)   =  g_o1d			! o1d
-        gind(i_o2)    =  g_o2			! o2
-        gind(i_o3)    =  g_o3			! o3
-        gind(i_h)     =  g_h			! h
-        gind(i_h2)    =  g_h2			! h2
-        gind(i_oh)    =  g_oh			! oh
-        gind(i_ho2)   =  g_ho2			! ho2
-        gind(i_h2o2)  =  g_h2o2			! h2o2
-        gind(i_n2)    =  g_n2			! n2
-        gind(i_ar)    =  g_ar			! ar
-        gind(i_h2o)   =  g_h2o_vap		! h2o
+
+         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
-c
-        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_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_n2)    = 1.034e3
-        cpi(i_ar)    = 1.000e3    !?
-        cpi(i_h2o)   = 1.870e3
-c
+c    Thermal conductivity and specific heat coefficients
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
-        nlayer=nlayermx
-        ngrid=ngridmx
-        nq=nqmx
+         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
 
-        firstcall=.false.
-      endif ! of if (firstcall)
+         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
 
-      DO l=1,nlayer
-        DO ig=1,ngrid
-          DO n=1,ncomptot
-            q2(ig,l,n)=max(1.e-30,
-     .         pq(ig,l,gind(n))+pdq(ig,l,gind(n))*ptimestep)
-          ENDDO
-          zt(ig,l)=pt(ig,l) +pdt(ig,l)*ptimestep
-        ENDDO
-      ENDDO
-  
-      do l=1,nlayermx
-        do ig=1,ngridmx
-          ntot=pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
-          cpnew(ig,l)=0.
-          akknew(ig,l)=0.
-          mmean(ig,l)=0.
-          q2tot(ig,l)=0.
-          nt=0.
-          do n=1,ncomptot
-            ni(n)=0.0 
-            do k=1,ncomptot
-              if(k.ne.n) ni(n)=ni(n)+q2(ig,l,k)/mmol(gind(k))
-            enddo
-            ni(n)=ntot/(1.+mmol(gind(n))/q2(ig,l,n)*ni(n))
-            cpnew(ig,l)=cpnew(ig,l)+ni(n)*cpi(n)
-            akknew(ig,l)=akknew(ig,l)+ni(n)*aki(n)
-            mmean(ig,l)=mmean(ig,l)+q2(ig,l,n)/mmol(gind(n))
-            q2tot(ig,l)=q2tot(ig,l)+q2(ig,l,n)
-c        if(ig.eq.1.and.l.eq.1) write(*,*)'q2tot(ig,l)',n,q2tot(ig,l) 
-            if(cpi(n) .ne. 0.0) nt=nt+ni(n)
-          enddo 
+         firstcall=.false.
 
-c     print*,"concentrations rep 3",l,ig,nt,mmean(ig,l),zt(ig+1,l)
+      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
 
-          cpnew(ig,l)=cpnew(ig,l)/nt
-          akknew(ig,l)=akknew(ig,l)/nt
-          mmean(ig,l)=1/mmean(ig,l)    ! in amu
-          rnew(ig,l)=8.314/mmean(ig,l)*1.e3     ! J/kg/K		
-        enddo
-c        print*,l,mmean(1,l),cpnew(1,l),rnew(1,l)
-c         write(228,*),l,pplay(1,l),ntot
-      enddo
       return
       end 
Index: trunk/LMDZ.MARS/libf/phymars/callradite.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 358)
+++ trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 370)
@@ -3,5 +3,4 @@
      $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
      &     tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice)
-     &     tauref,tau,aerosol,ccn,rdust,rice,nuice,co2ice)
 
        IMPLICIT NONE
