Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 606)
+++ 	(revision )
@@ -1,502 +1,0 @@
-      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
-     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
-     $                   dqscloud,tauref,co2ice,
-     $                   pu,pdu,pv,pdv,surfdust,surfice)
-c
-      implicit none
-c
-c=======================================================================
-c
-c   subject:
-c   --------
-c
-c  Prepare the call for the photochemical module, and send back the
-c  tendencies from photochemistry in the chemical species mass mixing ratios
-c
-c   Author:   Sebastien Lebonnois (08/11/2002)
-c   -------
-c    update 12/06/2003 for water ice clouds and compatibility with dust
-c    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
-c    update 03/05/2005 cosmetic changes (Franck Lefevre)
-c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
-c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
-c    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
-c
-c   Arguments:
-c   ----------
-c
-c  Input:
-c
-c    ptimestep                  timestep (s)
-c    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
-c    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
-c    pt(ngridmx,nlayermx)       Temperature (K)
-c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
-c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
-c    pdu(ngridmx,nlayermx)      u component tendency (K)
-c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
-c    pdv(ngridmx,nlayermx)      v component tendency (K)
-c    dist_sol                   distance of the sun (AU)
-c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
-c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
-c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
-c    tauref(ngridmx)            Optical depth at 7 hPa
-c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
-c    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
-c    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
-c
-c  Output:
-c
-c    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
-c    dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
-c
-c=======================================================================
-
-c    Declarations :
-c    --------------
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "chimiedata.h"
-#include "tracer.h"
-#include "comcstfi.h"
-#include "callkeys.h"
-#include "conc.h"
-
-c Arguments :
-c -----------
-
-c   inputs:
-c   -------
-
-      real    ptimestep
-      real    pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
-      real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
-      real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
-      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
-      real    pt(ngridmx,nlayermx)       ! temperature
-      real    pdt(ngridmx,nlayermx)      ! temperature tendency
-      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
-      real    pdu(ngridmx,nlayermx)      ! u component tendency
-      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
-      real    pdv(ngridmx,nlayermx)      ! v component tendency
-      real    dist_sol                   ! distance of the sun (AU)
-      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
-      real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
-      real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
-      real    zday                       ! date (time since Ls=0, in martian days)
-      real    tauref(ngridmx)            ! optical depth at 7 hPa
-      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
-      real    surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
-      real    surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
-
-c   outputs:
-c   --------
-
-      real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
-      real dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
-      real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
-      real dqscloud(ngridmx,nqmx)        ! tendencies on qsurf 
-
-c Local variables :
-c -----------------
-
-      character*20 str20
-      integer  ig,l,i,iq
-      integer  foundswitch, lswitch
-      integer  ig_vl1
-      real latvl1, lonvl1
-      real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
-                                     ! new mole fraction after
-      real zt(ngridmx,nlayermx)      ! temperature
-      real zu(ngridmx,nlayermx)      ! u component of the wind
-      real zv(ngridmx,nlayermx)      ! v component of the wind
-      real taucol                    ! optical depth at 7 hPa
-
-      logical depos                  ! switch for dry deposition
-      parameter (depos=.false.)
-c
-c  for each column of atmosphere:
-c
-      real zpress(nlayermx)       !  Pressure (mbar)
-      real zdens(nlayermx)        !  Density  (cm-3)
-      real ztemp(nlayermx)        !  Temperature (K)
-      real zlocal(nlayermx)       !  Altitude (km)
-      real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
-      real szacol                 !  Solar zenith angle
-      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
-      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
-      real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
-c
-c for output:
-c
-      real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
-
-      logical output              ! to issue calls to writediagfi and stats
-      parameter (output=.true.)   ! see at end of routine
-
-      logical,save :: firstcall=.true.
-      integer,save :: nbq       ! number of tracers used in the chemistry
-      integer,save :: niq(nqmx) ! array storing the indexes of the tracers
-
-! index of tracers:
-
-      integer,save :: i_co2=0
-      integer,save :: i_co=0
-      integer,save :: i_o=0
-      integer,save :: i_o1d=0
-      integer,save :: i_o2=0
-      integer,save :: i_o3=0
-      integer,save :: i_h=0
-      integer,save :: i_h2=0
-      integer,save :: i_oh=0
-      integer,save :: i_ho2=0
-      integer,save :: i_h2o2=0
-      integer,save :: i_ch4=0
-      integer,save :: i_n2=0
-      integer,save :: i_ar=0
-      integer,save :: i_ice=0 ! water ice
-      integer,save :: i_h2o=0 ! water vapour
-c
-c=======================================================================
-c     initialization of the chemistry (first call only)
-c=======================================================================
-c
-      if (firstcall) then
-c
-         if (photochem) then
-            print*,'calchim: Read photolysis lookup table'
-            call read_phototable
-         end if
-
-         ! find index of chemical tracers to use
-         nbq=0 ! to count number of tracers
-         i_co2=igcm_co2
-         if (i_co2.eq.0) then
-           write(*,*) "calchim: Error; no CO2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_co2
-         endif
-         i_co=igcm_co
-         if (i_co.eq.0) then
-           write(*,*) "calchim: Error; no CO tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_co
-         endif
-         i_o=igcm_o
-         if (i_o.eq.0) then
-           write(*,*) "calchim: Error; no O tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_o
-         endif
-         i_o1d=igcm_o1d
-         if (i_o1d.eq.0) then
-           write(*,*) "calchim: Error; no O1D tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_o1d
-         endif
-         i_o2=igcm_o2
-         if (i_o2.eq.0) then
-           write(*,*) "calchim: Error; no O2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_o2
-         endif
-         i_o3=igcm_o3
-         if (i_o3.eq.0) then
-           write(*,*) "calchim: Error; no O3 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_o3
-         endif
-         i_h=igcm_h
-         if (i_h.eq.0) then
-           write(*,*) "calchim: Error; no H tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_h
-         endif
-         i_h2=igcm_h2
-         if (i_h2.eq.0) then
-           write(*,*) "calchim: Error; no H2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_h2
-         endif
-         i_oh=igcm_oh
-         if (i_oh.eq.0) then
-           write(*,*) "calchim: Error; no OH tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_oh
-         endif
-         i_ho2=igcm_ho2
-         if (i_ho2.eq.0) then
-           write(*,*) "calchim: Error; no HO2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_ho2
-         endif
-         i_h2o2=igcm_h2o2
-         if (i_h2o2.eq.0) then
-           write(*,*) "calchim: Error; no H2O2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_h2o2
-         endif
-         i_ch4=igcm_ch4
-         if (i_ch4.eq.0) then
-           write(*,*) "calchim: Error; no CH4 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_ch4
-         endif
-         i_n2=igcm_n2
-         if (i_n2.eq.0) then
-           write(*,*) "calchim: Error; no N2 tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_n2
-         endif
-         i_ar=igcm_ar
-         if (i_ar.eq.0) then
-           write(*,*) "calchim: Error; no AR tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_ar
-         endif
-         i_ice=igcm_h2o_ice
-         if (i_ice.eq.0) then
-           write(*,*) "calchim: Error; no water ice tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_ice
-         endif
-         i_h2o=igcm_h2o_vap
-         if (i_h2o.eq.0) then
-           write(*,*) "calchim: Error; no water vapor tracer !!!"
-           stop
-         else
-           nbq=nbq+1
-           niq(nbq)=i_h2o
-         endif
-
-         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
-         write(*,*) '               i_co2  = ',i_co2
-         write(*,*) '               i_co   = ',i_co
-         write(*,*) '               i_o    = ',i_o
-         write(*,*) '               i_o1d  = ',i_o1d
-         write(*,*) '               i_o2   = ',i_o2
-         write(*,*) '               i_o3   = ',i_o3
-         write(*,*) '               i_h    = ',i_h
-         write(*,*) '               i_h2   = ',i_h2
-         write(*,*) '               i_oh   = ',i_oh
-         write(*,*) '               i_ho2  = ',i_ho2
-         write(*,*) '               i_h2o2 = ',i_h2o2
-         write(*,*) '               i_ch4  = ',i_ch4
-         write(*,*) '               i_n2   = ',i_n2
-         write(*,*) '               i_ar   = ',i_ar
-         write(*,*) '               i_ice  = ',i_ice
-         write(*,*) '               i_h2o  = ',i_h2o
-         
-         firstcall = .false.
-      end if ! if (firstcall)
-
-! Initialize output tendencies to zero (to handle case of tracers which
-! are not used in the chemistry (e.g. dust))
-
-      dqchim(:,:,:) = 0
-      dqschim(:,:)  = 0
-c
-!     latvl1= 22.27
-!     lonvl1= -47.94
-!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
-!    $        int(1.5+(lonvl1+180)*iim/360.)
-c
-c=======================================================================
-c     loop over grid
-c=======================================================================
-c
-      do ig = 1,ngridmx
-c
-c     local updates
-c 
-         foundswitch = 0
-         do l = 1,nlayermx
-            do i = 1,nbq
-              iq = niq(i) ! get tracer index
-              zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-              zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-            end do
-            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
-            zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
-            zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
-c
-            zpress(l) = pplay(ig,l)/100.
-            ztemp(l)  = zt(ig,l)
-            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
-            zlocal(l) = zzlay(ig,l)/1000.
-c
-c           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
-c
-            surfdust1d(l) = surfdust(ig,l)*1.e-2
-            surfice1d(l)  = surfice(ig,l)*1.e-2
-c
-c           search for switch index between regions
-c
-            if (photochem .and. thermochem) then
-               if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
-                  lswitch = l
-                  foundswitch=1
-               end if
-            end if
-            if ( .not. photochem) then
-               lswitch = 22
-            end if
-            if (.not.  thermochem) then
-               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
-            end if
-c
-         end do ! of do l=1,nlayermx
-c
-         szacol = acos(mu0(ig))*180./pi
-         taucol = tauref(ig)
-c
-c=======================================================================
-c     call chemical subroutine
-c=======================================================================
-c
-c        chemistry in lower atmosphere
-c
-         if (photochem) then
-            call photochemistry(lswitch,zycol,szacol,ptimestep,
-     $                          zpress,ztemp,zdens,dist_sol,
-     $                          surfdust1d,surfice1d,jo3,taucol)
-         end if
-c
-c        chemistry in upper atmosphere
-c
-         if (thermochem) then
-            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
-     $                       zlocal,szacol,ptimestep,zday)
-         end if
-c
-c        dry deposition
-c
-         if (depos) then
-            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev, 
-     $                      zu, zv, zt, zycol, ptimestep, co2ice)
-         end if
-c
-c=======================================================================
-c     tendencies
-c=======================================================================
-c
-c     must be 0. for water ice:
-c
-         if (water) then
-            do l = 1,nlayermx
-               dqchim(ig,l,i_ice) = 0.
-            end do
-         end if
-c
-c     tendency for CO2 = - sum of others for lower atmosphere
-c     tendency for O   = - sum of others for upper atmosphere
-c
-         do l = 1,nlayermx
-            if (l .lt. lswitch) then
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then
-                     dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)
-     $                                - zq(ig,l,iq))/ptimestep
-                  else if (iq.eq.i_co2) then
-                     dqchim(ig,l,iq) = 0.
-                  end if
-                  dqschim(ig,iq) = 0.
-               end do ! of do i=1,nbq
-               
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  if (iq.ne.i_co2) then
-                     dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2) 
-     $                                  - dqchim(ig,l,iq)
-                  end if
-               end do
-             else if (l .ge. lswitch) then
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                   if ((iq.ne.i_o).and.(iq.ne.i_ice)) then
-                      dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)
-     $                                  /mmean(ig,l)
-     $                                 - zq(ig,l,iq))/ptimestep
-                   else if (iq.eq.i_o) then
-                      dqchim(ig,l,iq) = 0.
-                   end if
-               enddo
-               
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                   if (iq.ne.i_o) then
-                      dqchim(ig,l,i_o) = dqchim(ig,l,i_o) 
-     $                                 - dqchim(ig,l,iq)
-                   end if
-               end do
-             end if ! of if (l.lt.lswitch) else if (l.ge.lswitch)
-          end do ! of do l = 1,nlayermx
-c
-c     condensation of h2o2
-c
-          call perosat(ig,ptimestep,pplev,pplay,
-     $                 ztemp,zycol,dqcloud,dqscloud)
-c
-c     density and j(o3->o1d), for outputs
-c
-          do l = 1,nlayermx
-             jo3_3d(ig,l) = jo3(l)
-          end do
-c
-c=======================================================================
-c     end of loop over grid
-c=======================================================================
-c
-      end do ! of do ig=1,ngridmx
-c
-c=======================================================================
-c     write outputs
-c=======================================================================
-c
-! value of parameter 'output' to trigger writting of outputs
-! is set above at the declaration of the variable.
-
-      if (output) then
-         if (ngridmx .gt. 1) then
-            call writediagfi(ngridmx,'jo3','j o3->o1d',
-     $                       's-1',3,jo3_3d(1,1))
-           if (callstats) then
-            call wstats(ngridmx,'jo3','j o3->o1d',
-     $                       's-1',3,jo3_3d(1,1))
-           endif
-         end if ! of if (ngridmx.gt.1)
-      endif ! of if (output)
-c
-      end
Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 618)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 618)
@@ -0,0 +1,420 @@
+      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
+                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
+                         dqscloud,tauref,co2ice,                            &
+                         pu,pdu,pv,pdv,surfdust,surfice)
+
+      implicit none
+
+!=======================================================================
+!
+!   subject:
+!   --------
+!
+!  Prepare the call for the photochemical module, and send back the
+!  tendencies from photochemistry in the chemical species mass mixing ratios
+!
+!   Author:   Sebastien Lebonnois (08/11/2002)
+!   -------
+!    update 12/06/2003 for water ice clouds and compatibility with dust
+!    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
+!    update 03/05/2005 cosmetic changes (Franck Lefevre)
+!    update sept. 2008 identify tracers by their names (Ehouarn Millour)
+!    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
+!    update 16/03/2012 optimization (Franck Lefevre)
+!
+!   Arguments:
+!   ----------
+!
+!  Input:
+!
+!    ptimestep                  timestep (s)
+!    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
+!    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
+!    pt(ngridmx,nlayermx)       Temperature (K)
+!    pdt(ngridmx,nlayermx)      Temperature tendency (K)
+!    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
+!    pdu(ngridmx,nlayermx)      u component tendency (K)
+!    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
+!    pdv(ngridmx,nlayermx)      v component tendency (K)
+!    dist_sol                   distance of the sun (AU)
+!    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
+!    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
+!    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
+!    tauref(ngridmx)            Optical depth at 7 hPa
+!    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
+!    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
+!    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
+!
+!  Output:
+!
+!    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
+!    dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
+!
+!=======================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "chimiedata.h"
+#include "tracer.h"
+#include "comcstfi.h"
+#include "callkeys.h"
+#include "conc.h"
+
+!     input:
+
+      real :: ptimestep
+      real :: pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
+      real :: zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
+      real :: pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
+      real :: zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
+      real :: pt(ngridmx,nlayermx)       ! temperature
+      real :: pdt(ngridmx,nlayermx)      ! temperature tendency
+      real :: pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
+      real :: pdu(ngridmx,nlayermx)      ! u component tendency
+      real :: pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
+      real :: pdv(ngridmx,nlayermx)      ! v component tendency
+      real :: dist_sol                   ! distance of the sun (AU)
+      real :: mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
+      real :: pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
+      real :: pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
+      real :: zday                       ! date (time since Ls=0, in martian days)
+      real :: tauref(ngridmx)            ! optical depth at 7 hPa
+      real :: co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
+      real :: surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
+      real :: surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
+
+!     output:
+
+      real :: dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
+      real :: dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
+      real :: dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
+      real :: dqscloud(ngridmx,nqmx)        ! tendencies on qsurf 
+
+!     local variables:
+
+      integer,save :: nbq                   ! number of tracers used in the chemistry
+      integer,save :: niq(nqmx)             ! array storing the indexes of the tracers
+      integer :: major(nlayermx)            ! index of major species
+      integer :: ig,l,i,iq,iqmax
+      integer :: foundswitch, lswitch
+
+      integer,save :: i_co2  = 0
+      integer,save :: i_co   = 0
+      integer,save :: i_o    = 0
+      integer,save :: i_o1d  = 0
+      integer,save :: i_o2   = 0
+      integer,save :: i_o3   = 0
+      integer,save :: i_h    = 0
+      integer,save :: i_h2   = 0
+      integer,save :: i_oh   = 0
+      integer,save :: i_ho2  = 0
+      integer,save :: i_h2o2 = 0
+      integer,save :: i_ch4  = 0
+      integer,save :: i_n2   = 0
+      integer,save :: i_h2o  = 0
+
+      integer :: ig_vl1
+
+      real    :: latvl1, lonvl1
+      real    :: zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
+                                           ! new mole fraction after
+      real    :: zt(ngridmx,nlayermx)      ! temperature
+      real    :: zu(ngridmx,nlayermx)      ! u component of the wind
+      real    :: zv(ngridmx,nlayermx)      ! v component of the wind
+      real    :: taucol                    ! optical depth at 7 hPa
+
+      logical,save :: firstcall = .true.
+      logical,save :: depos = .false.      ! switch for dry deposition
+
+!     for each column of atmosphere:
+
+      real :: zpress(nlayermx)       !  Pressure (mbar)
+      real :: zdens(nlayermx)        !  Density  (cm-3)
+      real :: ztemp(nlayermx)        !  Temperature (K)
+      real :: zlocal(nlayermx)       !  Altitude (km)
+      real :: zycol(nlayermx,nqmx)   !  Composition (mole fractions)
+      real :: szacol                 !  Solar zenith angle
+      real :: surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
+      real :: surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
+      real :: jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
+
+!     for output:
+
+      logical :: output                 ! to issue calls to writediagfi and stats
+      parameter (output = .true.)
+      real :: jo3_3d(ngridmx,nlayermx)  ! Photodissociation rate O3->O1D (s-1)
+
+!=======================================================================
+!     initialization of the chemistry (first call only)
+!=======================================================================
+
+      if (firstcall) then
+
+         if (photochem) then
+            print*,'calchim: Read photolysis lookup table'
+            call read_phototable
+         end if
+
+         ! find index of chemical tracers to use
+         nbq = 0 ! to count number of tracers
+         i_co2 = igcm_co2
+         if (i_co2 == 0) then
+            write(*,*) "calchim: Error; no CO2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_co2
+         end if
+         i_co = igcm_co
+         if (i_co == 0) then
+            write(*,*) "calchim: Error; no CO tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_co
+         end if
+         i_o = igcm_o
+         if (i_o == 0) then
+            write(*,*) "calchim: Error; no O tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_o
+         end if
+         i_o1d = igcm_o1d
+         if (i_o1d == 0) then
+            write(*,*) "calchim: Error; no O1D tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_o1d
+         end if
+         i_o2 = igcm_o2
+         if (i_o2 == 0) then
+            write(*,*) "calchim: Error; no O2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_o2
+         end if
+         i_o3 = igcm_o3
+         if (i_o3 == 0) then
+            write(*,*) "calchim: Error; no O3 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_o3
+         end if
+         i_h = igcm_h
+         if (i_h == 0) then
+            write(*,*) "calchim: Error; no H tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_h
+         end if
+         i_h2 = igcm_h2
+         if (i_h2 == 0) then
+            write(*,*) "calchim: Error; no H2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_h2
+         end if
+         i_oh = igcm_oh
+         if (i_oh == 0) then
+            write(*,*) "calchim: Error; no OH tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_oh
+         end if
+         i_ho2 = igcm_ho2
+         if (i_ho2 == 0) then
+            write(*,*) "calchim: Error; no HO2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_ho2
+         end if
+         i_h2o2 = igcm_h2o2
+         if (i_h2o2 == 0) then
+            write(*,*) "calchim: Error; no H2O2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_h2o2
+         end if
+         i_ch4 = igcm_ch4
+         if (i_ch4 == 0) then
+            write(*,*) "calchim: no CH4 tracer !!!"
+            write(*,*) "CH4 will be ignored in the chemistry"
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_ch4
+         end if
+         i_n2 = igcm_n2
+         if (i_n2 == 0) then
+            write(*,*) "calchim: Error; no N2 tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_n2
+         end if
+         i_h2o = igcm_h2o_vap
+         if (i_h2o == 0) then
+            write(*,*) "calchim: Error; no water vapor tracer !!!"
+            stop
+         else
+            nbq = nbq + 1
+            niq(nbq) = i_h2o
+         end if
+
+         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
+         
+         firstcall = .false.
+      end if ! if (firstcall)
+
+! Initializations
+
+      zycol(:,:)    = 0.
+      dqchim(:,:,:) = 0
+      dqschim(:,:)  = 0
+
+!     latvl1= 22.27
+!     lonvl1= -47.94
+!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +    &
+!             int(1.5+(lonvl1+180)*iim/360.)
+
+!=======================================================================
+!     loop over grid
+!=======================================================================
+
+      do ig = 1,ngridmx
+
+         foundswitch = 0
+         do l = 1,nlayermx
+            do i = 1,nbq
+               iq = niq(i) ! get tracer index
+               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
+               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+            end do
+            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
+            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
+            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
+            zpress(l) = pplay(ig,l)/100.
+            ztemp(l)  = zt(ig,l)
+            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
+            zlocal(l) = zzlay(ig,l)/1000.
+
+!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
+
+            surfdust1d(l) = surfdust(ig,l)*1.e-2
+            surfice1d(l)  = surfice(ig,l)*1.e-2
+
+!           search for switch index between regions
+
+            if (photochem .and. thermochem) then
+               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
+                  lswitch = l
+                  foundswitch = 1
+               end if
+            end if
+            if (.not. photochem) then
+               lswitch = 22
+            end if
+            if (.not. thermochem) then
+               lswitch = min(50,nlayermx+1)
+            end if
+
+         end do ! of do l=1,nlayermx
+
+         szacol = acos(mu0(ig))*180./pi
+         taucol = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
+
+!=======================================================================
+!     call chemical subroutines
+!=======================================================================
+
+!        chemistry in lower atmosphere
+
+         if (photochem) then
+            call photochemistry(lswitch,zycol,szacol,ptimestep,    &
+                                zpress,ztemp,zdens,dist_sol,       &
+                                surfdust1d,surfice1d,jo3,taucol)
+
+!        ozone photolysis, for output
+
+            do l = 1,nlayermx
+               jo3_3d(ig,l) = jo3(l)
+            end do
+
+!        condensation of h2o2
+
+            call perosat(ig,ptimestep,pplev,pplay,                 &
+                         ztemp,zycol,dqcloud,dqscloud)
+         end if
+
+!        chemistry in upper atmosphere
+
+         if (thermochem) then
+            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,  &
+                             zlocal,szacol,ptimestep,zday)
+         end if
+
+!        dry deposition
+
+         if (depos) then
+            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,& 
+                            zu, zv, zt, zycol, ptimestep, co2ice)
+         end if
+
+!=======================================================================
+!     tendencies
+!=======================================================================
+
+!     index of the most abundant species at each level
+
+         major(:) = maxloc(zycol, dim = 2)
+
+!     tendency for the most abundant species = - sum of others
+
+         do l = 1,nlayermx
+            iqmax = major(l)
+            do i = 1,nbq
+               iq = niq(i) ! get tracer index
+               if (iq /= iqmax) then
+                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)  &
+                                   - zq(ig,l,iq))/ptimestep
+                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax)              &
+                                     - dqchim(ig,l,iq) 
+               end if
+            end do
+         end do ! of do l = 1,nlayermx
+
+!=======================================================================
+!     end of loop over grid
+!=======================================================================
+
+      end do ! of do ig=1,ngridmx
+
+!=======================================================================
+!     write outputs
+!=======================================================================
+
+! value of parameter 'output' to trigger writting of outputs
+! is set above at the declaration of the variable.
+
+      if (photochem .and. output) then
+         if (ngridmx > 1) then
+            call writediagfi(ngridmx,'jo3','j o3->o1d',    &
+                             's-1',3,jo3_3d(1,1))
+           if (callstats) then
+              call wstats(ngridmx,'jo3','j o3->o1d',       &
+                          's-1',3,jo3_3d(1,1))
+           endif
+         end if ! of if (ngridmx.gt.1)
+      end if ! of if (output)
+
+      return
+      end
Index: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 606)
+++ trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 618)
@@ -3,17 +3,16 @@
       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=======================================================================
+!=======================================================================
+! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
+!
+! mmean(ngridmx,nlayermx)	amu
+! cpnew(ngridmx,nlayermx)	J/kg/K
+! rnew(ngridmx,nlayermx)	J/kg/K
+! akknew(ngridmx,nlayermx)	coefficient of thermal concduction
+!
+! version: April 2012 - Franck Lefevre
+!=======================================================================
 
-c    Declarations
-c    ------------
+!     declarations
  
 #include "dimensions.h"
@@ -26,6 +25,5 @@
 #include "conc.h"
 
-c    Input/Output
-c    ------------
+!     input/output
 
       real pplay(ngridmx,nlayermx)
@@ -36,109 +34,109 @@
       real ptimestep
 
-c    Local variables
-c    ---------------
+!     local variables
 
-      integer       :: l, ig, n, k
-      integer, save :: gind(ncomp)
+      integer       :: i, l, ig, iq
+      integer, save :: nbq, niq(nqmx)
       real          :: ni(nqmx), ntot
-      real          :: zq(ngridmx,nlayermx,ncomp)
-      real          :: zt(ngridmx,nlayermx)
-      real, save    :: aki(ncomp)
-      real, save    :: cpi(ncomp)
+      real          :: zq(ngridmx, nlayermx, nqmx)
+      real          :: zt(ngridmx, nlayermx)
+      real, save    :: aki(nqmx)
+      real, save    :: cpi(nqmx)
 
       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
+!        find index of chemical tracers to use
+!        initialize thermal conductivity and specific heat coefficients
+!        !? values are estimated
 
-         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
+         nbq = 0 ! to count number of tracers used in this subroutine
 
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c    Thermal conductivity and specific heat coefficients
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+         if (igcm_co2 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_co2
+            aki(nbq) = 3.072e-4
+            cpi(nbq) = 0.834e3
+         end if
+         if (igcm_co /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_co
+            aki(nbq) = 4.87e-4
+            cpi(nbq) = 1.034e3
+         end if
+         if (igcm_o /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_o
+            aki(nbq) = 7.59e-4
+            cpi(nbq) = 1.3e3
+         end if
+         if (igcm_o1d /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_o1d
+            aki(nbq) = 7.59e-4  !?
+            cpi(nbq) = 1.3e3    !?
+         end if
+         if (igcm_o2 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_o2
+            aki(nbq) = 5.68e-4
+            cpi(nbq) = 0.9194e3
+         end if
+         if (igcm_o3 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_o3
+            aki(nbq) = 3.00e-4  !?
+            cpi(nbq) = 0.800e3  !?
+         end if
+         if (igcm_h /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_h
+            aki(nbq) = 0.0
+            cpi(nbq) = 20.780e3
+         end if
+         if (igcm_h2 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_h2
+            aki(nbq) = 36.314e-4
+            cpi(nbq) = 14.266e3
+         end if
+         if (igcm_oh /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_oh
+            aki(nbq)  = 7.00e-4 !?
+            cpi(nbq)  = 1.045e3
+         end if
+         if (igcm_ho2 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_ho2
+            aki(nbq) = 0.0
+            cpi(nbq) = 1.065e3  !?
+         end if
+         if (igcm_n2 /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_n2
+            aki(nbq) = 5.6e-4
+            cpi(nbq) = 1.034e3
+         end if
+         if (igcm_ar /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_ar
+            aki(nbq) = 0.0      !?
+            cpi(nbq) = 1.000e3  !?
+         end if
+         if (igcm_h2o_vap /= 0) then
+            nbq = nbq + 1
+            niq(nbq) = igcm_h2o_vap
+            aki(nbq) = 0.0
+            cpi(nbq) = 1.870e3
+         end if
 
-         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.
 
-         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
+      end if ! if (firstcall)
 
-         firstcall=.false.
+!     update temperature
 
-      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
@@ -146,23 +144,27 @@
          end do
       end do
-c
-c     update tracers
-c
+
+!     update tracers
+
       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)
+            do i = 1,nbq
+               iq = niq(i) 
+               zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)
+     $                                 + pdq(ig,l,iq)*ptimestep)
             end do
          end do
       end do
-c
-c     mmean : mean molecular mass
-c     rnew  : specific gas constant
-c
+
+!     mmean : mean molecular mass
+!     rnew  : specific gas constant
+
+      mmean(:,:)  = 0.
+
       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))
+            do i = 1,nbq
+               iq = niq(i) 
+               mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
             end do
             mmean(ig,l) = 1./mmean(ig,l)
@@ -170,20 +172,24 @@
          end do
       end do
-c
-c     cpnew  : specicic heat
-c     akknew : thermal conductivity cofficient
-c      
+
+!     cpnew  : specicic heat
+!     akknew : thermal conductivity cofficient
+      
+      cpnew(:,:)  = 0.
+      akknew(:,:) = 0.
+
       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)
+            do i = 1,nbq
+               iq = niq(i) 
+               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
-c        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
+!        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
       end do
 
Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F	(revision 606)
+++ 	(revision )
@@ -1,644 +1,0 @@
-      subroutine inichim_newstart(pq,ps,sig,nq,latfi,lonfi,airefi,
-     $ thermo,qsurf)
-
-      implicit none
-
-c=======================================================================
-c
-c   subject:
-c   --------
-c
-c  Initialization of the composition for use in the building of a new start file
-c  used by program newstart.F
-c  also used by program testphys1d.F
-c
-c   Author:   Sebastien Lebonnois (08/11/2002)
-c   -------
-c   Revised 07/2003 by Monica Angelats-i-Coll to use input files
-c   Modified 10/2008 Identify tracers by their names Ehouarn Millour
-c   Modified 11/2011 Addition of methane (Franck Lefevre)
-c
-c   Arguments:
-c   ----------
-c
-c    pq(iip1,jjp1,llm,nqmx)  Advected fields, ie chemical species here
-c    sig = sigma                vertical coordinate (interface of layers)
-c    nq: number of tracers to initialize (used to evaluate if only
-c        chemistry species are to be initialized or if water vapour
-c        and water ice should also be initialized. 
-c        NB: number of chemistry tracers is ncomp (set in chimiedata.h)
-c=======================================================================
-
-c    Declarations :
-c    --------------
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "paramet.h"
-#include "chimiedata.h"
-#include "tracer.h"
-#include "comcstfi.h"
-#include "comdiurn.h"
-#include "callkeys.h"
-#include "temps.h"
-#include "datafile.h"
-
-c Arguments :
-c -----------
-
-      real    ps(iip1,jjp1)
-      real    pq(iip1,jjp1,llm,nqmx)     ! Advected fields, ie chemical species
-      real    sig(llm+1)                 ! vertical coordinate (interface of layers)
-      integer nq            !  =nqmx 
-                            ! or =nqmx-1 if H2O kept
-                            ! or =nqmx-2 if H2O and ice kept
-      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
-      integer thermo ! flag for thermosphere initialisation only
-      real :: qsurf(ngridmx,nqmx) ! surface tracers
-
-c Local variables :
-c -----------------
-
-      integer  iq,i,j,l,n, nspecini,inil,iqmax,iiq
-      INTEGER aa(nqmx)
-      REAL qtot(iip1,jjp1,llm)
-      real densconc,zgcm,zfile(252),vmrint,nt, mmean
-      real nxf(252),zfilemin(252),sigfile(252)
-      real vmrf(252,14),nf(252,14)
-      real tfile(252),zzfile(252)
-      real ni(14),niold(14)
-      real mmolf(14)    ! molecular mass in amu (species in files)
-      data mmolf/44.,40.,28.,32.,28.,16.,2.,   ! majors
-     &           1.,17.,33.,18.,34.,16.,48./   ! minors
-      character*3 tmp ! temporary variable
-      integer ierr
-
-      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
-      character(len=20) :: txt ! to store some text
-      integer :: nqchem_start,count
-      integer :: nqchem(nqmx) ! local chemistry tracer index array
-      integer :: nbqchem ! total number of chemical species (water included)
-
-c-----------------------------------------------------------------------
-c    Input/Output
-c    ------------
-! read in 'callphys.def'
-       call inichim_readcallphys()
-
-      ! check if tracers have 'old' names
-      oldnames=.false.
-      count=0
-      do iq=1,nqmx
-        txt=" "
-        write(txt,'(a1,i2.2)') 'q',iq
-        if (txt.eq.noms(iq)) then
-          count=count+1
-        endif
-      enddo ! of do iq=1,nqmx
-
-      if (count.eq.nqmx) then
-        write(*,*) "inichim_newstart: tracers seem to follow old ",
-     &             "naming convention (q01,q02,...)"
-        write(*,*) "   => will work for now ... "
-        write(*,*) "      but you should rename them"
-        oldnames=.true.
-      endif
-
-      ! copy/set tracer names
-      if (oldnames) then ! old names (derived from iq & option values)
-        ! 1. dust:
-        if (dustbin.ne.0) then ! transported dust
-          do iq=1,dustbin
-            txt=" "
-            write(txt,'(a4,i2.2)') 'dust',iq
-            noms(iq)=txt
-            mmol(iq)=100.
-          enddo
-          ! special case if doubleq
-          if (doubleq) then
-            if (dustbin.ne.2) then
-              write(*,*) 'initracer: error expected dustbin=2'
-            else
-!              noms(1)='dust_mass'   ! dust mass mixing ratio
-!              noms(2)='dust_number' ! dust number mixing ratio
-! same as above, but avoid explict possible out-of-bounds on array
-               noms(1)='dust_mass'   ! dust mass mixing ratio
-               do iq=2,2
-               noms(iq)='dust_number' ! dust number mixing ratio
-               enddo
-            endif
-          endif
-        endif
-        ! 2. water & ice
-        if (water) then
-          noms(nqmx)='h2o_vap'
-          mmol(nqmx)=18.
-!            noms(nqmx-1)='h2o_ice'
-!            mmol(nqmx-1)=18.
-          !"loop" to avoid potential out-of-bounds in array
-          do iq=nqmx-1,nqmx-1
-            noms(iq)='h2o_ice'
-            mmol(iq)=18.
-          enddo
-        endif
-        ! 3. Chemistry
-        if (photochem .or. callthermos) then
-          if (doubleq) then
-            nqchem_start=3
-          else
-            nqchem_start=dustbin+1
-          end if
-        endif ! of if (photochem .or. callthermos)
-        do iq=nqchem_start,nqchem_start+ncomp-2
-          noms(iq)=nomchem(iq-nqchem_start+1)
-          mmol(iq)=mmolchem(iq-nqchem_start+1)
-        enddo
-        ! 4. Other tracers ????
-        if ((dustbin.eq.0).and.(.not.water)) then
-          noms(1)='co2'
-          mmol(1)=44
-          if (nqmx.eq.2) then
-            noms(nqmx)='Ar_N2'
-            mmol(nqmx)=30
-          endif
-        endif
-      else
-        ! noms(:) already contain tracer names
-        ! mmol(:) array is initialized later (see below)
-      endif ! of if (oldnames)
-
-      ! special modification when using 'old' tracers:
-      ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour
-      ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice
-      if (oldnames.and.water) then
-        write(*,*)'inichim_newstart: moving surface water ice to index '
-     &            ,nqmx-1
-        ! "loop" to avoid potential out-of-bounds on arrays
-        do iq=nqmx-1,nqmx-1
-        qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
-        enddo
-        qsurf(1:ngridmx,nqmx)=0
-      endif 
-
-c Dimension test:
-c ---------------
-
-       if (water) then
-!         if ((nqchem_min+ncomp+1).ne.nqmx) then
-         if ((dustbin+ncomp+2).ne.nqmx) then
-!            print*,'********* Dimension problem! ********'
-!            print*,"nqchem_min+ncomp+1).ne.nqmx"
-            print*,'inichim_newstart: tracer dimension problem:'
-            print*,"(dustbin+ncomp+2).ne.nqmx"
-            print*,"ncomp: ",ncomp
-!            print*,"nqchem_min: ",nqchem_min
-            print*,"nqmx: ",nqmx
-            print*,'Change ncomp in chimiedata.h'
-            STOP
-         endif
-       else
-!         if ((nqchem_min+ncomp).ne.nqmx) then
-          if ((dustbin+ncomp+1).ne.nqmx) then
-!            print*,'********* Dimension problem! ********'
-!            print*,"nqchem_min+ncomp).ne.nqmx"
-            print*,'initracer: tracer dimension problem:'
-            print*,"(dustbin+ncomp+1).ne.nqmx"
-            print*,"ncomp: ",ncomp
-!            print*,"nqchem_min: ",nqchem_min
-            print*,"nqmx: ",nqmx
-            print*,'Change ncomp in chimiedata.h'
-            STOP
-         endif
-       endif
-
-c noms and mmol vectors:
-c ----------------------
-!
-!       if (iceparty) then
-!          do iq=nqchem_min,nqmx-2
-!            noms(iq) = nomchem(iq-nqchem_min+1)
-!            mmol(iq) = mmolchem(iq-nqchem_min+1)
-!          enddo
-!          noms(nqmx-1) = 'ice'
-!          mmol(nqmx-1) = 18.
-!          noms(nqmx)   = 'h2o'
-!          mmol(nqmx)   = 18.
-!       else
-!          do iq=nqchem_min,nqmx-1
-!            noms(iq) = nomchem(iq-nqchem_min+1)
-!            mmol(iq) = mmolchem(iq-nqchem_min+1)
-!          enddo
-!          noms(nqmx)   = 'h2o'
-!          mmol(nqmx)   = 18.
-!       endif
-!       if (nqchem_min.gt.1) then
-!          do iq=1,nqchem_min-1
-!            noms(iq) = 'dust'
-!            mmol(iq) = 100.
-!          enddo
-!       endif
-
-
-! Identify tracers by their names: (and set corresponding values of mmol)
-      ! 0. initialize tracer indexes to zero:
-      do iq=1,nqmx
-        igcm_dustbin(iq)=0
-      enddo
-      igcm_dust_mass=0
-      igcm_dust_number=0
-      igcm_h2o_vap=0
-      igcm_h2o_ice=0
-      igcm_co2=0
-      igcm_co=0
-      igcm_o=0
-      igcm_o1d=0
-      igcm_o2=0
-      igcm_o3=0
-      igcm_h=0
-      igcm_h2=0
-      igcm_oh=0
-      igcm_ho2=0
-      igcm_h2o2=0
-      igcm_ch4=0
-      igcm_n2=0
-      igcm_ar=0
-      igcm_ar_n2=0
-
-      ! 1. find dust tracers
-      count=0
-      if (dustbin.gt.0) then
-        do iq=1,nqmx
-          txt=" "
-          write(txt,'(a4,i2.2)')'dust',count+1
-          if (noms(iq).eq.txt) then
-            count=count+1
-            igcm_dustbin(count)=iq
-            mmol(iq)=100.
-          endif
-        enddo !do iq=1,nqmx
-      endif ! of if (dustbin.gt.0)
-      if (doubleq) then
-        do iq=1,nqmx
-          if (noms(iq).eq."dust_mass") then
-            igcm_dust_mass=iq
-            count=count+1
-          endif
-          if (noms(iq).eq."dust_number") then
-            igcm_dust_number=iq
-            count=count+1
-          endif
-        enddo
-      endif ! of if (doubleq)
-      ! 2. find chemistry and water tracers
-      nbqchem=0
-      do iq=1,nqmx
-        if (noms(iq).eq."co2") then
-          igcm_co2=iq
-          mmol(igcm_co2)=44.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."co") then
-          igcm_co=iq
-          mmol(igcm_co)=28.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."o") then
-          igcm_o=iq
-          mmol(igcm_o)=16.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."o1d") then
-          igcm_o1d=iq
-          mmol(igcm_o1d)=16.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."o2") then
-          igcm_o2=iq
-          mmol(igcm_o2)=32.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."o3") then
-          igcm_o3=iq
-          mmol(igcm_o3)=48.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."h") then
-          igcm_h=iq
-          mmol(igcm_h)=1.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."h2") then
-          igcm_h2=iq
-          mmol(igcm_h2)=2.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."oh") then
-          igcm_oh=iq
-          mmol(igcm_oh)=17.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."ho2") then
-          igcm_ho2=iq
-          mmol(igcm_ho2)=33.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."h2o2") then
-          igcm_h2o2=iq
-          mmol(igcm_h2o2)=34.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."ch4") then
-          igcm_ch4=iq
-          mmol(igcm_ch4)=16.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."n2") then
-          igcm_n2=iq
-          mmol(igcm_n2)=28.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."ar") then
-          igcm_ar=iq
-          mmol(igcm_ar)=40.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."h2o_vap") then
-          igcm_h2o_vap=iq
-          mmol(igcm_h2o_vap)=18.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        if (noms(iq).eq."h2o_ice") then
-          igcm_h2o_ice=iq
-          mmol(igcm_h2o_ice)=18.
-          count=count+1
-          nbqchem=nbqchem+1
-        endif
-        ! Other stuff: e.g. for simulations using co2 + neutral gaz
-        if (noms(iq).eq."Ar_N2") then
-          igcm_ar_n2=iq
-          mmol(igcm_ar_n2)=30.
-          count=count+1
-        endif
-      enddo ! of do iq=1,nqmx
-      
-      ! check that we identified all tracers:
-      if (count.ne.nqmx) then
-        write(*,*) "inichim_newstart: found only ",count," tracers"
-        write(*,*) "               expected ",nqmx
-        do iq=1,count
-          write(*,*)'      ',iq,' ',trim(noms(iq))
-        enddo
-        stop
-      else
-       write(*,*)"inichim_newstart: found all expected tracers, namely:"
-        do iq=1,nqmx
-          write(*,*)'      ',iq,' ',trim(noms(iq))
-        enddo
-      endif
-
-      ! build local chemistry tracer index array nqchem(:)
-      ! as in the old days, water vapour is last and water ice,
-      ! if present, is just before water vapour
-      do iq=1,16 ! loop so as to avoid out-of-bounds on array
-      if (iq==1) nqchem(i)=igcm_co2
-      if (iq==2) nqchem(i)=igcm_co
-      if (iq==3) nqchem(i)=igcm_o
-      if (iq==4) nqchem(i)=igcm_o1d
-      if (iq==5) nqchem(i)=igcm_o2
-      if (iq==6) nqchem(i)=igcm_o3
-      if (iq==7) nqchem(i)=igcm_h
-      if (iq==8) nqchem(i)=igcm_h2
-      if (iq==9) nqchem(i)=igcm_oh
-      if (iq==10) nqchem(i)=igcm_ho2
-      if (iq==11) nqchem(i)=igcm_h2o2
-      if (iq==12) nqchem(i)=igcm_ch4
-      if (iq==13) nqchem(i)=igcm_n2
-      if (iq==14) nqchem(i)=igcm_ar
-      if (iq==15) nqchem(i)=igcm_h2o_ice
-      if (iq==16) nqchem(i)=igcm_h2o_vap
-      enddo
-
-! Load in chemistry data for initialization:
-c----------------------------------------------------------------------
-c Order of Major species in file:
-c
-c    1=CO2 
-c    2=AR
-c    3=N2  
-c    4=O2  
-c    5=CO  
-c    6=O   
-c    7=H2
-c
-c Order of Minor species in files:
-c
-c    1=H  
-c    2=OH 
-c    3=HO2 
-c    4=H2O 
-c    5=H2O2 
-c    6=O1D
-c    7=O3
-c
-c----------------------------------------------------------------------
-c Major species:
-        aa(igcm_co2) = 1
-        aa(igcm_ar)  = 2
-        aa(igcm_n2)  = 3
-        aa(igcm_o2)  = 4
-        aa(igcm_co)  = 5
-        aa(igcm_o)   = 6
-        aa(igcm_h2)  = 7
-c Minor species:
-        aa(igcm_h)       = 8
-        aa(igcm_oh)      = 9
-        aa(igcm_ho2)     = 10
-        aa(igcm_h2o_vap) = 11
-        aa(igcm_h2o2)    = 12
-        aa(igcm_o1d)     = 13 
-        aa(igcm_o3)      = 14
-c----------------------------------------------------------------------
-      open(210, iostat=ierr,file=
-     & trim(datafile)//'/atmosfera_LMD_may.dat')
-      if (ierr.ne.0) then
-        write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
-        write(*,*)'(in aeronomars/inichim.F)'
-        write(*,*)'It should be in :', trim(datafile),'/'
-        write(*,*)'1) You can change this path in callphys.def with'
-        write(*,*)'   datadir=/path/to/datafiles/'
-        write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
-        write(*,*)'   can be obtained online on:'
-        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
-         STOP
-      endif
-      open(220, iostat=ierr,file=
-     & trim(datafile)//'/atmosfera_LMD_min.dat')
-      if (ierr.ne.0) then
-        write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
-        write(*,*)'(in aeronomars/inichim.F)'
-        write(*,*)'It should be in :', trim(datafile),'/'
-        write(*,*)'1) You can change this path in callphys.def with'
-        write(*,*)'   datadir=/path/to/datafiles/'
-        write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
-        write(*,*)'   can be obtained online on:'
-        write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
-         STOP
-      endif
-
-      read(210,*) tmp
-      read(220,*) tmp
-        do l=1,252
-          read(210,112) tfile(l),zfile(l),nxf(l),(vmrf(l,n),n=1,7)
-          zfile(l)=zfile(l)*100.              !pressure (Pa)
-          sigfile(l)=zfile(l)/zfile(1)       
-        enddo
-        do l=1,252
-          read(220,113)zfilemin(l),(vmrf(l,n),n=8,14)
-          zfilemin(l)=zfilemin(l)*1000.       !height (m)
-          do n=1,14
-            nf(l,n)=vmrf(l,n)*nxf(l)
-          enddo
-        enddo
-        close(210)
-        close(220)
-
-c flag thermo set to 1 or 0 in newstart 
-c inil=33 for initialisation of thermosphere only        
-c inil=1 for initialisation of all layers        
-        if (thermo.eq.1) then
-        inil=33
-        else
-        inil=1
-        endif
-
-! Initialization of tracers
-
-      do i=1,iip1
-       do j=1,jjp1
-        do l=inil,llm
-
-c          zgcm=sig(l)
-          zgcm=sig(l)*ps(i,j)
-          densconc=0.
-          nt=0.
-
-          do n=1,14
-            call intrplf(zgcm,vmrint,zfile,nf(1,n),252)
-c            call intrplf(zgcm,vmrint,sigfile,nf(1,n),252)
-            ni(n)=vmrint
-
-            densconc=densconc+ni(n)*mmolf(n)
-            nt=nt+ni(n)
-          enddo
-
-          mmean=densconc/nt                     ! in amu
-
-          if (nq.ne.nbqchem) then ! don't initialize water vapour
-            do n=1,nq-1
-            pq(i,j,l,nqchem(n))=
-     &                 ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean
-            if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n)) 
-            enddo
-            if (water .and. (nq.gt.nbqchem-2)) then
-              pq(i,j,l,nqchem(nq)) = 0.
-              if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) 
-            else
-              pq(i,j,l,nqchem(nq))=
-     &              ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
-              if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) 
-            endif
-          endif ! of if (nq.ne.nbqchem)
-          
-          if (nq.eq.nbqchem) then ! also initialize water vapour
-            do n=1,nq-2
-              pq(i,j,l,nqchem(n))=
-     &                   ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean
-              if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n)) 
-            enddo
-            if (water) then
-              pq(i,j,l,nqchem(nq-1)) = 0.
-              if(i.eq.iip1) then
-                pq(i,j,l,nqchem(nq-1))=pq(1,j,l,nqchem(nq-1)) 
-              endif
-              pq(i,j,l,nqchem(nq))=
-     &                   ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
-              if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) 
-            else
-              do n=nq-1,nq
-                pq(i,j,l,nqchem(nq))=
-     &                   ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
-                if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq)) 
-              enddo
-            endif ! of if (water)
-          endif ! of if (nq.eq.nqmx)
-        enddo     !nlayer
-       enddo      !ngrid_j
-      enddo       !ngrid_i
-
-cccccccccccccccccccccccccccccc      
-c  make sure that sum of q = 1      
-c dominent species is = 1 - sum(all other species)      
-cccccccccccccccccccccccccccccc      
-!      iqmax=nqchem_min
-       iqmax=1
-     
-!      if ((nqmx-nqchem_min).gt.10) then
-      if (nbqchem.gt.10) then
-       do l=1,llm
-        do j=1,jjp1
-          do i=1,iip1
-!           do iq=nqchem_min,nqmx
-            do iiq=1,nbqchem
-              iq=nqchem(iiq)
-              if ( (pq(i,j,l,iq).gt.pq(i,j,l,iqmax)).and.
-     .           (noms(iq).ne."ice") )  then 
-              iqmax=iq
-              endif
-            enddo
-            pq(i,j,l,iqmax)=1.
-            qtot(i,j,l)=0
-!           do iq=nqchem_min,nqmx
-            do iiq=1,nbqchem
-             iq=nqchem(iiq)
-             if ( (iq.ne.iqmax).and.
-     .           (noms(iq).ne."ice") ) then        
-               pq(i,j,l,iqmax)=pq(i,j,l,iqmax)-pq(i,j,l,iq)        
-             endif
-            enddo !iq
-!            do iq=nqchem_min,nqmx
-            do iiq=1,nbqchem
-             iq=nqchem(iiq)
-              if (noms(iq).ne."ice") then
-                qtot(i,j,l)=qtot(i,j,l)+pq(i,j,l,iq)
-              endif
-c            if (i.eq.1.and.j.eq.1.and.l.eq.1) write(*,*) 'qtot(i,j,l)',
-c     $    qtot(i,j,l)
-            enddo !iq
-            if (i.eq.1.and.j.eq.1) write(*,*) 'inichim_newstart: ',
-     $        'qtot(i,j,l)=',qtot(i,j,l)
-          enddo !i   
-         enddo !j   
-       enddo !l  
-      endif ! of if (nbqchem.gt.10)
-ccccccccccccccccccccccccccccccc
-
-66      format(i2,6(1x,e11.5))
-112     format(7x,f7.3,3x,e12.6,3x,e12.6,7(3x,e12.6))
-113     format(1x,f6.2,7(3x,e12.6))
-
-      end
Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 618)
+++ trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 618)
@@ -0,0 +1,530 @@
+      subroutine inichim_newstart(pq, qsurf, ps, flagh2o, flagthermo)
+
+      implicit none
+
+!=======================================================================
+!
+!  Purpose:
+!  --------
+!
+!  Initialization of the chemistry for use in the building of a new start file
+!  used by program newstart.F
+!  also used by program testphys1d.F
+!
+!  Authors: 
+!  -------
+!  Initial version 11/2002 by Sebastien Lebonnois
+!  Revised 07/2003 by Monica Angelats-i-Coll to use input files
+!  Modified 10/2008 Identify tracers by their names Ehouarn Millour
+!  Modified 11/2011 Addition of methane Franck Lefevre
+!  Rewritten 04/2012 Franck Lefevre
+!
+!  Arguments:
+!  ----------
+!
+!  pq(iip1,jjp1,llm,nqmx)  Advected fields, ie chemical species here
+!  qsurf(ngridmx,nqmx)     Amount of tracer on the surface (kg/m2)
+!  ps(iip1,jjp1)           Surface pressure (Pa)
+!  flagh2o                 flag for initialisation of h2o (1: yes / 0: no)
+!  flagthermo              flag for initialisation of thermosphere only (1: yes / 0: no)
+!
+!=======================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "paramet.h"
+#include "tracer.h"
+#include "comvert.h"
+#include "callkeys.h"
+#include "datafile.h"
+
+! inputs :
+
+      real,intent(in) :: ps(iip1,jjp1)            ! surface pressure in the gcm (Pa)   
+      integer,intent(in) :: flagh2o               ! flag for h2o initialisation
+      integer,intent(in) :: flagthermo            ! flag for thermosphere initialisation only
+
+! outputs :
+
+      real,intent(out) :: pq(iip1,jjp1,llm,nqmx)  ! advected fields, ie chemical species
+      real,intent(out) :: qsurf(ngridmx,nqmx)     ! surface values (kg/m2) of tracers
+
+! local :
+
+      integer :: iq, i, j, l, n, nbqchem
+      integer :: count, ierr, dummy
+      real    :: mmean(iip1,jjp1,llm)             ! mean molecular mass (g)
+      real    :: pgcm                             ! pressure at each layer in the gcm (Pa)
+
+      integer, parameter         :: nalt = 252    ! number of altitudes in the initialization files
+      integer, parameter         :: nspe = 14     ! number of species in the initialization files
+      integer, dimension(nspe)   :: niq           ! local index of species in initialization files
+      real, dimension(nalt)      :: tinit, zzfile ! temperature in initialization files
+      real, dimension(nalt)      :: pinit         ! pressure in initialization files
+      real, dimension(nalt)      :: densinit      ! total number density in initialization files
+      real, dimension(nalt,nspe) :: vmrinit       ! mixing ratios in initialization files
+      real, dimension(nspe)      :: vmrint        ! mixing ratio interpolated onto the gcm vertical grid
+      real                       :: vmr
+
+      character(len=20)          :: txt           ! to store some text
+
+! 1. identify tracers by their names: (and set corresponding values of mmol)
+
+! 1.1 initialize tracer indexes to zero:
+
+      do iq = 1,nqmx
+        igcm_dustbin(iq) = 0
+      end do
+
+      igcm_dust_mass   = 0
+      igcm_dust_number = 0
+      igcm_ccn_mass    = 0
+      igcm_ccn_number  = 0
+      igcm_h2o_vap     = 0
+      igcm_h2o_ice     = 0
+      igcm_co2         = 0
+      igcm_co          = 0
+      igcm_o           = 0
+      igcm_o1d         = 0
+      igcm_o2          = 0
+      igcm_o3          = 0
+      igcm_h           = 0
+      igcm_h2          = 0
+      igcm_oh          = 0
+      igcm_ho2         = 0
+      igcm_h2o2        = 0
+      igcm_ch4         = 0
+      igcm_n2          = 0
+      igcm_ar          = 0
+      igcm_ar_n2       = 0
+      igcm_co2plus     = 0
+      igcm_oplus       = 0
+      igcm_o2plus      = 0
+      igcm_coplus      = 0
+      igcm_cplus       = 0
+      igcm_nplus       = 0
+      igcm_noplus      = 0
+      igcm_n2plus      = 0
+      igcm_hplus       = 0
+      igcm_elec        = 0
+
+! 1.2 find dust tracers
+
+      count = 0
+
+      if (dustbin > 0) then
+         do iq = 1,nqmx
+            txt = " "
+            write(txt,'(a4,i2.2)') 'dust', count + 1
+            if (noms(iq) == txt) then
+               count = count + 1
+               igcm_dustbin(count) = iq
+               mmol(iq) = 100.
+            end if
+         end do !do iq=1,nqmx
+      end if ! of if (dustbin.gt.0)
+
+      if (doubleq) then
+         do iq = 1,nqmx
+            if (noms(iq) == "dust_mass") then
+               igcm_dust_mass = iq
+               count = count + 1
+            end if
+            if (noms(iq) == "dust_number") then
+               igcm_dust_number = iq
+               count = count + 1
+            end if
+         end do
+      end if ! of if (doubleq)
+
+      if (scavenging) then
+         do iq = 1,nqmx
+            if (noms(iq) == "ccn_mass") then
+               igcm_ccn_mass = iq
+               count = count + 1
+            end if
+            if (noms(iq) == "ccn_number") then
+               igcm_ccn_number = iq
+               count = count + 1
+            end if
+         end do
+      end if ! of if (scavenging)
+
+      if (submicron) then
+         do iq=1,nqmx
+            if (noms(iq) == "dust_submicron") then
+               igcm_dust_submicron = iq
+               mmol(iq) = 100.
+               count = count + 1
+            end if
+         end do
+      end if ! of if (submicron)
+
+! 1.3 find chemistry and water tracers
+
+      nbqchem = 0
+      do iq = 1,nqmx
+         if (noms(iq) == "co2") then
+            igcm_co2 = iq
+            mmol(igcm_co2) = 44.
+            count = count + 1
+            nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "co") then
+           igcm_co = iq
+           mmol(igcm_co) = 28.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "o") then
+           igcm_o = iq
+           mmol(igcm_o) = 16.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "o1d") then
+           igcm_o1d = iq
+           mmol(igcm_o1d) = 16.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "o2") then
+           igcm_o2 = iq
+           mmol(igcm_o2) = 32.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "o3") then
+           igcm_o3 = iq
+           mmol(igcm_o3) = 48.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "h") then
+           igcm_h = iq
+           mmol(igcm_h) = 1.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "h2") then
+           igcm_h2 = iq
+           mmol(igcm_h2) = 2.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "oh") then
+           igcm_oh = iq
+           mmol(igcm_oh) = 17.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "ho2") then
+           igcm_ho2 = iq
+           mmol(igcm_ho2) = 33.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "h2o2") then
+           igcm_h2o2 = iq
+           mmol(igcm_h2o2) = 34.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "ch4") then
+           igcm_ch4 = iq
+           mmol(igcm_ch4) = 16.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "n2") then
+           igcm_n2 = iq
+           mmol(igcm_n2) = 28.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "n") then
+           igcm_n = iq
+           mmol(igcm_n) = 14.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "n2d") then
+           igcm_n2d = iq
+           mmol(igcm_n2d) = 14.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "no") then
+           igcm_no = iq
+           mmol(igcm_no) = 30.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "no2") then
+           igcm_no2 = iq
+           mmol(igcm_no2) = 46.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "ar") then
+           igcm_ar = iq
+           mmol(igcm_ar) = 40.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "h2o_vap") then
+           igcm_h2o_vap = iq
+           mmol(igcm_h2o_vap) = 18.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "h2o_ice") then
+           igcm_h2o_ice = iq
+           mmol(igcm_h2o_ice) = 18.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+
+! 1.4 find ions
+
+        if (noms(iq) == "co2plus") then
+           igcm_co2plus = iq
+           mmol(igcm_co2plus) = 44.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "oplus") then
+           igcm_oplus = iq
+           mmol(igcm_oplus) = 16.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "o2plus") then
+           igcm_o2plus = iq
+           mmol(igcm_o2plus) = 32.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "coplus") then
+           igcm_coplus = iq
+           mmol(igcm_coplus) = 28.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "cplus") then
+           igcm_cplus = iq
+           mmol(igcm_cplus) = 12.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "nplus") then
+           igcm_nplus = iq
+           mmol(igcm_nplus) = 14.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "noplus") then
+           igcm_noplus = iq
+           mmol(igcm_noplus) = 30.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "n2plus") then
+           igcm_n2plus = iq
+           mmol(igcm_n2plus) = 28.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "hplus") then
+           igcm_hplus = iq
+           mmol(igcm_hplus) = 1.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
+        if (noms(iq) == "elec") then
+           igcm_elec = iq
+           mmol(igcm_elec) = 1./1822.89
+           count = count + 1
+        end if
+
+! 1.5 find idealized non-condensible tracer
+
+        if (noms(iq) == "Ar_N2") then
+           igcm_ar_n2 = iq
+           mmol(igcm_ar_n2) = 30.
+           count = count + 1
+        end if
+
+      end do ! of do iq=1,nqmx
+      
+! 1.6 check that we identified all tracers:
+
+      if (count /= nqmx) then
+         write(*,*) "inichim_newstart: found only ",count," tracers"
+         write(*,*) "                  expected ",nqmx
+         do iq = 1,count
+            write(*,*) '      ', iq, ' ', trim(noms(iq))
+         end do
+         stop
+      else
+         write(*,*) "inichim_newstart: found all expected tracers"
+         do iq = 1,nqmx
+            write(*,*) '      ', iq, ' ', trim(noms(iq))
+         end do
+      end if
+
+! 2. load in chemistry data for initialization:
+
+! order of major species in initialization file:
+!
+!    1: co2 
+!    2: ar
+!    3: n2  
+!    4: o2  
+!    5: co  
+!    6: o   
+!    7: h2
+!
+! order of minor species in initialization file:
+!
+!    1: h  
+!    2: oh 
+!    3: ho2 
+!    4: h2o 
+!    5: h2o2 
+!    6: o1d
+!    7: o3
+
+! major species:
+
+      niq(1) = igcm_co2
+      niq(2) = igcm_ar
+      niq(3) = igcm_n2
+      niq(4) = igcm_o2
+      niq(5) = igcm_co
+      niq(6) = igcm_o
+      niq(7) = igcm_h2
+
+! minor species:
+
+      niq(8)  = igcm_h
+      niq(9)  = igcm_oh
+      niq(10) = igcm_ho2
+      niq(11) = igcm_h2o_vap
+      niq(12) = igcm_h2o2
+      niq(13) = igcm_o1d 
+      niq(14) = igcm_o3
+
+! 2.1 open initialization files
+
+      open(210, iostat=ierr,file=trim(datafile)// &
+                            '/atmosfera_LMD_may.dat')
+      if (ierr /= 0) then
+         write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
+         write(*,*)'(in aeronomars/inichim.F)'
+         write(*,*)'It should be in :', trim(datafile),'/'
+         write(*,*)'1) You can change this path in callphys.def with'
+         write(*,*)'   datadir=/path/to/datafiles/'
+         write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
+         write(*,*)'   can be obtained online on:'
+         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
+         stop
+      end if
+      open(220, iostat=ierr,file=trim(datafile)// &
+                            '/atmosfera_LMD_min.dat')
+      if (ierr /= 0) then
+         write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
+         write(*,*)'(in aeronomars/inichim.F)'
+         write(*,*)'It should be in :', trim(datafile),'/'
+         write(*,*)'1) You can change this path in callphys.def with'
+         write(*,*)'   datadir=/path/to/datafiles/'
+         write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
+         write(*,*)'   can be obtained online on:'
+         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
+         stop
+      end if
+
+! 2.2 read initialization files
+
+! major species
+
+      read(210,*)
+      do l = 1,nalt
+         read(210,*) dummy, tinit(l), pinit(l), densinit(l), &
+                     (vmrinit(l,n), n = 1,7)
+         pinit(l) = pinit(l)*100.              ! conversion in Pa
+         pinit(l) = log(pinit(l))              ! for the vertical interpolation
+      end do
+      close(210)
+
+! minor species
+
+      read(220,*)
+      do l = 1,nalt
+         read(220,*) dummy, (vmrinit(l,n), n = 8,14)
+      end do 
+      close(220)
+
+! 3. initialization of tracers
+
+      do i = 1,iip1
+         do j = 1,jjp1
+            do l = 1,llm
+
+               pgcm = aps(l) + bps(l)*ps(i,j)  ! gcm pressure
+               pgcm = log(pgcm)                ! for the vertical interpolation
+               mmean(i,j,l) = 0.
+
+! 3.1 vertical interpolation
+
+               do n = 1,nspe
+                  call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt)
+                  vmrint(n) = vmr
+                  iq = niq(n)
+                  mmean(i,j,l) = mmean(i,j,l) + vmrint(n)*mmol(iq)
+               end do
+
+! 3.2 attribute mixing ratio: - all layers or only thermosphere
+!                             - with our without h2o 
+
+               if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 1.e-3)) then
+                  do n = 1,nspe
+                     iq = niq(n)
+                     if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
+                        pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
+                     end if
+                  end do
+               end if
+
+            end do
+         end do
+      end do
+
+      ! set surface values of chemistry tracers to zero
+      if (flagthermo == 0) then
+        ! NB: no problem for "surface water vapour" tracer which is always 0
+        do n=1,nspe
+          iq=niq(n)
+          qsurf(1:ngridmx,iq)=0
+        enddo
+      endif
+
+
+! 3.3 initialization of tracers not contained in the initialization files
+
+! methane : 10 ppbv
+
+      if (igcm_ch4 /= 0) then
+         vmr = 10.e-9       
+         do i = 1,iip1
+            do j = 1,jjp1
+               do l = 1,llm
+                  pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l)
+               end do
+            end do
+         end do
+         ! set surface value to zero
+         qsurf(1:ngridmx,igcm_ch4)=0
+      end if
+
+      end
Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_readcallphys.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/inichim_readcallphys.F	(revision 606)
+++ 	(revision )
@@ -1,520 +1,0 @@
-      SUBROUTINE inichim_readcallphys
-!
-!=======================================================================
-!
-!   subject:
-!   --------
-!
-!   Initialisation for the physical parametrisations of the LMD 
-!   martian atmospheric general circulation model.
-!   This routine is called by inichim_newstart (which is used by programs
-!   newstart and testphys1d but not the GCM)
-!
-!   author: Frederic Hourdin 15 / 10 /93
-!   -------
-!   modifications: Sebastien Lebonnois 11/06/2003 (new callphys.def)
-!                  10/2008 adapted to using tracer sby name. Ehouarn Millour
-!                  07/2009 use 'getin' to read callphys.def. EM
-!
-!   arguments:
-!   ----------
-!
-!   input:
-!   ------
-!
-!=======================================================================
-!
-!-----------------------------------------------------------------------
-!   declarations:
-!   -------------
-! to use  'getin'
-      USE ioipsl_getincom 
-      IMPLICIT NONE
-#include "dimensions.h"
-#include "dimphys.h"
-#include "planete.h"
-#include "comcstfi.h"
-#include "comsaison.h"
-#include "comdiurn.h"
-#include "comgeomfi.h"
-#include "callkeys.h"
-#include "surfdat.h"
-
-      character*12 ch1
-      integer ierr
-      logical chem, h2o
-
-
-! --------------------------------------------------------------
-!  Reading the "callphys.def" file controlling some key options
-! --------------------------------------------------------------
-     
-      ! check that 'callphys.def' file is around
-      OPEN(99,file='callphys.def',status='old',form='formatted'
-     &     ,iostat=ierr)
-      CLOSE(99)
-      
-      IF(ierr.EQ.0) THEN
-         PRINT*
-         PRINT*
-         PRINT*,'--------------------------------------------'
-         PRINT*,' inichim_readcallphys: Parametres pour la physique',
-     &          ' (callphys.def)'
-         PRINT*,'--------------------------------------------'
-
-
-         write(*,*) "Run with or without tracer transport ?"
-         tracer=.false. ! default value
-         call getin("tracer",tracer)
-         write(*,*) " tracer = ",tracer
-
-         write(*,*) "Diurnal cycle ?"
-         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
-         diurnal=.true. ! default value
-         call getin("diurnal",diurnal)
-         write(*,*) " diurnal = ",diurnal
-
-         write(*,*) "Seasonal cycle ?"
-         write(*,*) "(if season=False, Ls stays constant, to value ",
-     &   "set in 'start'"
-         season=.true. ! default value
-         call getin("season",season)
-         write(*,*) " season = ",season
-
-         write(*,*) "Write some extra output to the screen ?"
-         lwrite=.false. ! default value
-         call getin("lwrite",lwrite)
-         write(*,*) " lwrite = ",lwrite
-
-         write(*,*) "Save statistics in file stats.nc ?"
-         callstats=.true. ! default value
-         call getin("callstats",callstats)
-         write(*,*) " callstats = ",callstats
-
-         write(*,*) "Save EOF profiles in file 'profiles' for ",
-     &              "Climate Database?"
-         calleofdump=.false. ! default value
-         call getin("calleofdump",calleofdump)
-         write(*,*) " calleofdump = ",calleofdump
-
-         write(*,*) "Dust scenario:"
-         iaervar=3 ! default value
-         call getin("iaervar",iaervar)
-         write(*,*) " iaervar = ",iaervar
-
-         write(*,*) "Dust vertical distribution:"
-         write(*,*) "(=1 Dust opt.deph read in startfi;",
-     & " =2 Viking scenario; =3 MGS scenario,",
-     & " =4 Mars Year 24 from TES assimilation)"
-         iddist=3 ! default value
-         call getin("iddist",iddist)
-         write(*,*) " iddist = ",iddist
-
-         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
-         topdustref= 90.0 ! default value
-         call getin("topdustref",topdustref)
-         write(*,*) " topdustref = ",topdustref
-
-
-         write(*,*) "call radiative transfer ?"
-         callrad=.true. ! default value
-         call getin("callrad",callrad)
-         write(*,*) " callrad = ",callrad
-
-         write(*,*) "call NLTE radiative schemes ?",
-     &              "(matters only if callrad=T)"
-         callnlte=.false. ! default value
-         call getin("callnlte",callnlte)
-         write(*,*) " callnlte = ",callnlte
-         
-         write(*,*) "call CO2 NIR absorption ?",
-     &              "(matters only if callrad=T)"
-         callnirco2=.false. ! default value
-         call getin("callnirco2",callnirco2)
-         write(*,*) " callnirco2 = ",callnirco2
-         
-         write(*,*) "call turbulent vertical diffusion ?"
-         calldifv=.true. ! default value
-         call getin("calldifv",calldifv)
-         write(*,*) " calldifv = ",calldifv
-
-         write(*,*) "call convective adjustment ?"
-         calladj=.true. ! default value
-         call getin("calladj",calladj)
-         write(*,*) " calladj = ",calladj
-         
-
-         write(*,*) "call CO2 condensation ?"
-         callcond=.true. ! default value
-         call getin("callcond",callcond)
-         write(*,*) " callcond = ",callcond
-
-         write(*,*)"call thermal conduction in the soil ?"
-         callsoil=.true. ! default value
-         call getin("callsoil",callsoil)
-         write(*,*) " callsoil = ",callsoil
-         
-
-         write(*,*)"call Lott's gravity wave/subgrid topography ",
-     &             "scheme ?"
-         calllott=.true. ! default value
-         call getin("calllott",calllott)
-         write(*,*)" calllott = ",calllott
-
-
-         write(*,*)"rad.transfer is computed every iradia",
-     &             " physical timestep"
-         iradia=1 ! default value
-         call getin("iradia",iradia)
-         write(*,*)" iradia = ",iradia
-         
-
-         write(*,*)"Output of the exchange coefficient mattrix ?",
-     &             "(for diagnostics only)"
-         callg2d=.false. ! default value
-         call getin("callg2d",callg2d)
-         write(*,*)" callg2d = ",callg2d
-
-         write(*,*)"Rayleigh scattering : (should be .false. for now)"
-         rayleigh=.false.
-         call getin("rayleigh",rayleigh)
-         write(*,*)" rayleigh = ",rayleigh
-
-
-! TRACERS:
-
-         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
-         dustbin=0 ! default value
-         call getin("dustbin",dustbin)
-         write(*,*)" dustbin = ",dustbin
-
-         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
-         active=.false. ! default value
-         call getin("active",active)
-         write(*,*)" active = ",active
-
-! Test of incompatibility:
-! if active is used, then dustbin should be > 0
-
-         if (active.and.(dustbin.lt.1)) then
-           print*,'if active is used, then dustbin should > 0'
-           stop
-         endif
-
-         write(*,*)"use mass and number mixing ratios to predict",
-     &             " dust size ?"
-         doubleq=.false. ! default value
-         call getin("doubleq",doubleq)
-         write(*,*)" doubleq = ",doubleq
-
-! Test of incompatibility:
-! if doubleq is used, then dustbin should be at least 2
-
-         if (doubleq.and.(dustbin.lt.2)) then
-           print*,'if doubleq is used, then dustbin should be > 1'
-           stop
-         endif
-
-         write(*,*)"dust lifted by GCM surface winds ?"
-         lifting=.false. ! default value
-         call getin("lifting",lifting)
-         write(*,*)" lifting = ",lifting
-
-! Test of incompatibility:
-! if lifting is used, then dustbin should be > 0
-
-         if (lifting.and.(dustbin.lt.1)) then
-           print*,'if lifting is used, then dustbin should > 0'
-           stop
-         endif
-
-         write(*,*)" dust lifted by dust devils ?"
-         callddevil=.false. !default value
-         call getin("callddevil",callddevil)
-         write(*,*)" callddevil = ",callddevil
-         
-
-! Test of incompatibility:
-! if dustdevil is used, then dustbin should be > 0
-
-         if (callddevil.and.(dustbin.lt.1)) then
-           print*,'if dustdevil is used, then dustbin should > 0'
-           stop
-         endif
-
-         write(*,*)"Dust scavenging by CO2 snowfall ?"
-         scavenging=.false. ! default value
-         call getin("scavenging",scavenging)
-         write(*,*)" scavenging = ",scavenging
-         
-
-! Test of incompatibility:
-! if scavenging is used, then dustbin should be > 0
-
-         if (scavenging.and.(dustbin.lt.1)) then
-           print*,'if scavenging is used, then dustbin should > 0'
-           stop
-         endif
-
-         write(*,*) "Gravitationnal sedimentation ?"
-         sedimentation=.true. ! default value
-         call getin("sedimentation",sedimentation)
-         write(*,*) " sedimentation = ",sedimentation
-
-         write(*,*) "includes water ice",
-     &              "(if true, 'water' must also be .true.)"
-
-         write(*,*) "Radiatively active transported atmospheric ",
-     &              "water ice ?"
-         activice=.false. ! default value
-         call getin("activice",activice)
-         write(*,*) " activice = ",activice
-
-         write(*,*) "Compute water cycle ?"
-         water=.false. ! default value
-         call getin("water",water)
-         write(*,*) " water = ",water
-         write(*,*) "Permanent water caps at poles ?",
-     &               " .true. is RECOMMENDED"
-         write(*,*) "(with .true., North cap is a source of water ",
-     &   "and South pole is a cold trap)"
-         caps=.true. ! default value
-         call getin("caps",caps)
-         write(*,*) " caps = ",caps
-
-! Test of incompatibility:
-! if activice is used, then water should be used too
-
-         if (activice.and..not.water) then
-           print*,'if activice is used, water should be used too'
-           stop
-         endif
-
-         write(*,*) "photochemistry: include chemical species"
-         photochem=.false. ! default value
-         call getin("photochem",photochem)
-         write(*,*) " photochem = ",photochem
-
-
-! THERMOSPHERE
-
-         write(*,*) "call thermosphere ?"
-         callthermos=.false. ! default value
-         call getin("callthermos",callthermos)
-         write(*,*) " callthermos = ",callthermos
-         
-
-         write(*,*) " water included without cycle ",
-     &              "(only if water=.false.)"
-         thermoswater=.false. ! default value
-         call getin("thermoswater",thermoswater)
-         write(*,*) " thermoswater = ",thermoswater
-
-         write(*,*) "call thermal conduction ?",
-     &    " (only if callthermos=.true.)"
-         callconduct=.false. ! default value
-         call getin("callconduct",callconduct)
-         write(*,*) " callconduct = ",callconduct
-
-         write(*,*) "call EUV heating ?",
-     &   " (only if callthermos=.true.)"
-         calleuv=.false.  ! default value
-         call getin("calleuv",calleuv)
-         write(*,*) " calleuv = ",calleuv
-
-         write(*,*) "call molecular viscosity ?",
-     &   " (only if callthermos=.true.)"
-         callmolvis=.false. ! default value
-         call getin("callmolvis",callmolvis)
-         write(*,*) " callmolvis = ",callmolvis
-
-         write(*,*) "call molecular diffusion ?",
-     &   " (only if callthermos=.true.)"
-         callmoldiff=.false. ! default value
-         call getin("callmoldiff",callmoldiff)
-         write(*,*) " callmoldiff = ",callmoldiff
-         
-
-         write(*,*) "call thermospheric photochemistry ?",
-     &   " (only if callthermos=.true.)"
-         thermochem=.false. ! default value
-         call getin("thermochem",thermochem)
-         write(*,*) " thermochem = ",thermochem
-
-         write(*,*) "date for solar flux calculation:",
-     &   " (1985 < date < 2002)"
-         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
-         solarcondate=1993.4 ! default value
-         call getin("solarcondate",solarcondate)
-         write(*,*) " solarcondate = ",solarcondate
-         
-
-         if (.not.callthermos) then
-           if (thermoswater) then
-             print*,'if thermoswater is set, callthermos must be true'
-             stop
-           endif          
-           if (callconduct) then
-             print*,'if callconduct is set, callthermos must be true'
-             stop
-           endif        
-           if (calleuv) then
-             print*,'if calleuv is set, callthermos must be true'
-             stop
-           endif         
-           if (callmolvis) then
-             print*,'if callmolvis is set, callthermos must be true'
-             stop
-           endif        
-           if (callmoldiff) then
-             print*,'if callmoldiff is set, callthermos must be true'
-             stop
-           endif          
-           if (thermochem) then
-             print*,'if thermochem is set, callthermos must be true'
-             stop
-           endif          
-        endif
-
-! Test of incompatibility:
-! if photochem is used, then water should be used too
-
-         if (photochem.and..not.water) then
-           print*,'if photochem is used, water should be used too'
-           stop
-         endif
-
-! if callthermos is used, then thermoswater should be used too 
-! (if water not used already)
-
-         if (callthermos .and. .not.water) then
-           if (callthermos .and. .not.thermoswater) then
-             print*,'if callthermos is used, water or thermoswater 
-     &               should be used too'
-             stop
-           endif
-         endif
-
-         PRINT*,'--------------------------------------------'
-         PRINT*
-         PRINT*
-      ELSE
-         write(*,*)
-         write(*,*) 'Cannot read file callphys.def. Is it here ?'
-         stop
-      ENDIF
-      CLOSE(99)
-
-      pi=2.*asin(1.)
-
-!     managing the tracers, and tests:
-!     -------------------------------
-
-      if(tracer) then
-
-!          when photochem is used, nqchem_min is the rank
-!          of the first chemical species
-
-! Ehouarn: nqchem_min is now meaningless and no longer used
-!       nqchem_min = 1
-       if (photochem .or. callthermos) then
-         chem = .true.
-!        if (doubleq) then
-!          nqchem_min = 3
-!        else
-!          nqchem_min = dustbin+1
-!        end if
-       end if
-
-       if (water .or. thermoswater) h2o = .true.
-
-c          TESTS
-
-       print*,'inichim_readcallphys: TRACERS:'
-
-       if ((doubleq).and.(h2o).and.
-     $     (chem)) then
-!         print*,' 1: dust ; 2: dust (doubleq)'
-!         print*,' 3 to ',nqmx-2,': chemistry'
-!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
-         print*,' 2 dust tracers (doubleq)'
-         print*,' 1 water vapour tracer'
-         print*,' 1 water ice tracer'
-         print*,nqmx-4,' chemistry tracers'
-       endif
-
-       if ((doubleq).and.(h2o).and.
-     $     .not.(chem)) then
-!         print*,' 1: dust ; 2: dust (doubleq)'
-!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
-         print*,' 2 dust tracers (doubleq)'
-         print*,' 1 water vapour tracer'
-         print*,' 1 water ice tracer'
-         if (nqmx.lt.4) then
-           print*,'nqmx should be at least 4 with these options.'
-		   print*,'(or check callphys.def)'
-           stop
-         endif
-       endif
-
-!       if ((doubleq).and..not.(h2o)) then
-!         print*,' 1: dust ; 2: dust (doubleq)'
-!         print*,' 2 dust tracers (doubleq)'
-!         if (nqmx.ne.2) then
-!           print*,'nqmx should be 2 with these options...'
-!		   print*,'(or check callphys.def)'
-!           stop
-!         endif
-!       endif
-
-!       if (.not.(doubleq).and.(h2o).and.
-!     $     (chem)) then
-!         if (dustbin.gt.0) then
-!           print*,' 1 to ',dustbin,': dust bins'
-!           print*,dustbin,' dust bins'
-!         endif
-!         print*,nqchem_min,' to ',nqmx-2,': chemistry'
-!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
-!         print*,nqmx-2-dustbin,' chemistry tracers'
-!         print*,' 1 water vapour tracer'
-!         print*,' 1 water ice tracer'
-!       endif
-
-!       if (.not.(doubleq).and.(h2o).and.
-!     $     .not.(chem)) then
-!         if (dustbin.gt.0) then
-!           print*,' 1 to ',dustbin,': dust bins'
-!           print*,dustbin,' dust bins'
-!         endif
-!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
-!         print*,' 1 water vapour tracer'
-!         print*,' 1 water ice tracer'
-!         if (nqmx.ne.(dustbin+2)) then
-!           print*,'nqmx should be ',(dustbin+2),
-!     $            ' with these options...'
-!		   print*,'(or check callphys.def)'
-!           stop
-!         endif
-!       endif
-
-!       if (.not.(doubleq).and..not.(h2o)) then
-!         if (dustbin.gt.0) then
-!           print*,' 1 to ',dustbin,': dust bins'
-!           print*,dustbin,' dust bins'
-!           if (nqmx.ne.dustbin) then
-!             print*,'nqmx should be ',dustbin,
-!     $              ' with these options...'
-!             print*,'(or check callphys.def)'
-!             stop
-!           endif
-!         else
-!           print*,'dustbin=',dustbin,
-!     $            ': tracer should be F with these options...'
-!     $           ,'UNLESS you just want to move tracers around '
-!         endif
-!       endif
-
-      endif ! of if (tracer)
-
-      RETURN
-      END
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 606)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 618)
@@ -60,5 +60,6 @@
       parameter (nesp = 16)      ! number of species in the chemistry
 c
-      real ctimestep             ! chemistry timestep (s) 
+      real stimestep             ! standard timestep for the chemistry (s)
+      real ctimestep             ! real timestep for the chemistry (s) 
       real rm(nlayermx,nesp)     ! species volume mixing ratio 
       real j(nlayermx,nd)        ! interpolated photolysis rates (s-1)
@@ -85,9 +86,12 @@
 c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     ctimestep : chemistry timestep (s)                         c
-c                 1/3 of physical timestep                       c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      phychemrat = 3
+c     stimestep  : standard timestep for the chemistry (s)       c
+c     ctimestep  : real timestep for the chemistry (s)           c
+c     phychemrat : ratio physical/chemical timestep              c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      stimestep = 600. ! standard value : 10 mn
+c
+      phychemrat = nint(ptimestep/stimestep)
 c
       ctimestep = ptimestep/real(phychemrat)
@@ -1175,5 +1179,4 @@
 c
       integer l,iq
-      logical,save :: firstcall = .true.
       
 c     tracer indexes in the chemistry:
@@ -1196,73 +1199,4 @@
       integer,parameter :: i_ox   = 16
 c      
-c     first call initializations
-c
-      if (firstcall) then
-
-c       identify the indexes of the tracers we need
-
-        if (igcm_co2.eq.0) then
-          write(*,*) "concentrations: Error; no CO2 tracer !!!"
-          stop
-        endif
-        if (igcm_co.eq.0) then
-          write(*,*) "concentrations: Error; no CO tracer !!!"
-          stop
-        endif
-        if (igcm_o.eq.0) then
-          write(*,*) "concentrations: Error; no O tracer !!!"
-          stop
-        endif
-        if (igcm_o1d.eq.0) then
-          write(*,*) "concentrations: Error; no O1D tracer !!!"
-          stop
-        endif
-        if (igcm_o2.eq.0) then
-          write(*,*) "concentrations: Error; no O2 tracer !!!"
-          stop
-        endif
-        if (igcm_o3.eq.0) then
-          write(*,*) "concentrations: Error; no O3 tracer !!!"
-          stop
-        endif
-        if (igcm_h.eq.0) then
-          write(*,*) "concentrations: Error; no H tracer !!!"
-          stop
-        endif
-        if (igcm_h2.eq.0) then
-          write(*,*) "concentrations: Error; no H2 tracer !!!"
-          stop
-        endif
-        if (igcm_oh.eq.0) then
-          write(*,*) "concentrations: Error; no OH tracer !!!"
-          stop
-        endif
-        if (igcm_ho2.eq.0) then
-          write(*,*) "concentrations: Error; no HO2 tracer !!!"
-          stop
-        endif
-        if (igcm_h2o2.eq.0) then
-          write(*,*) "concentrations: Error; no H2O2 tracer !!!"
-          stop
-        endif
-        if (igcm_ch4.eq.0) then
-          write(*,*) "concentrations: Error; no CH4 tracer !!!"
-          stop
-        endif
-        if (igcm_n2.eq.0) then
-          write(*,*) "concentrations: Error; no N2 tracer !!!"
-          stop
-        endif
-        if (igcm_ar.eq.0) then
-          write(*,*) "concentrations: Error; no AR tracer !!!"
-          stop
-        endif
-        if (igcm_h2o_vap.eq.0) then
-          write(*,*) "concentrations: Error; no water vapor tracer !!!"
-          stop
-        endif
-        firstcall = .false.
-      end if ! of if (firstcall)
-c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     initialise chemical species
@@ -1281,8 +1215,17 @@
          rm(l,i_ho2)  = max(zycol(l, igcm_ho2),     1.e-30)
          rm(l,i_h2o2) = max(zycol(l, igcm_h2o2),    1.e-30)
-         rm(l,i_ch4)  = max(zycol(l, igcm_ch4),     1.e-30)
          rm(l,i_n2)   = max(zycol(l, igcm_n2),      1.e-30)
          rm(l,i_h2o)  = max(zycol(l, igcm_h2o_vap), 1.e-30)
       end do 
+
+      if (igcm_ch4 .eq. 0) then
+         do l = 1,lswitch-1
+            rm(l,i_ch4) = 0.
+         end do
+      else
+         do l = 1,lswitch-1
+            rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
+         end do
+      end if
 c
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -1370,8 +1313,13 @@
          zycol(l, igcm_ho2)     = rm(l,i_ho2) 
          zycol(l, igcm_h2o2)    = rm(l,i_h2o2)
-         zycol(l, igcm_ch4)     = rm(l,i_ch4)
          zycol(l, igcm_n2)      = rm(l,i_n2)
          zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
       end do 
+
+      if (igcm_ch4 .ne. 0) then
+         do l = 1,lswitch-1
+            zycol(l,igcm_ch4) = rm(l,i_ch4)
+         end do
+      end if
 c
       return
@@ -2000,3 +1948,2 @@
       return
       end
-
