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
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,lnblnk
      external lnblnk

      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-1
          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_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."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,15 ! 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_n2
      if (iq==13) nqchem(i)=igcm_ar
      if (iq==14) nqchem(i)=igcm_h2o_ice
      if (iq==15) 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=
     & datafile(1:lnblnk(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 :', datafile(1:lnblnk(datafile)),'/'
        write(*,*)'1) You can change this directory address in '
        write(*,*)'   file phymars/datafile.h'
        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=
     & datafile(1:lnblnk(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 :', datafile(1:lnblnk(datafile)),'/'
        write(*,*)'1) You can change this directory address in '
        write(*,*)'   file phymars/datafile.h'
        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