subroutine read_mixratios(nmolec,nlev,P,x,molec_names)
      implicit none
c
c     Purpose: to read Venus' atmosphere "mixratios.VIRA" file
c
c     Outputs:
c       + nmolec: number of molecules
c       + nlev: number of atmospheric levels mixratios are defined for
c       + P: array of pressure levels [Pa]
c       + x(i,mol): molar concentration (mixratio) of species 'mol', at level 'i'
c       + molec_names: array containing the names of the nmolec molecules
c
      include 'max.inc'
      include 'formats.inc'
      integer nmolec,nlev,ios,nmolec_mf,i,mol,dim
      double precision x(1:Nmax,1:Nmol_max)
      double precision P(1:Nmax)
      character*(Nchar_mx) file_mixratios
      character*(Nchar_mx) strtmp1,strtmp2
      character*(Nchar_mx) string
      character*10 molec_names(1:Nmol_max)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine read_mixratios'

      file_mixratios='./mixratios.VIRA'
      open(11,file=file_mixratios(1:strlen(file_mixratios)),
     &        status='old',iostat=ios)
      if (ios.ne.0) then        ! file not found
         call error(label)
         write(*,*) 'Data file could not be found:'
         write(*,*) file_mixratios(1:strlen(file_mixratios))
         stop
      endif
      read(11,*)
      read(11,35) strtmp1,dim,strtmp2,nlev
      nmolec=dim-1
      read(11,*)
      read(11,'(a)') string
c     Debug
c      write(*,*) 'string="',string(1:strlen(string)),'"'
c      stop
c     Debug
      call identify_molecules(string,nmolec_mf,molec_names)
c     Debug
c      write(*,*) 'This is :',label(1:strlen(label))
c      do mol=1,nmolec_mf
c         write(*,*) 'molec_names(',mol,')=',molec_names(mol)
c      enddo
c      stop
c     Debug

      if (nmolec_mf.ne.nmolec) then
         call error(label)
         write(*,*) 'File: ',file_mixratios(1:strlen(file_mixratios)),
     &        ' is supposed to contain ',nmolec,' molecules'
         stop
      endif
      read(11,*)
      do i=1,nlev
         read(11,*) P(i),(x(i,mol),mol=1,nmolec)
c     convert pressures from hPa to Pa
         P(i)=P(i)*1.0D+2
      enddo
      close(11)

      return
      end



      subroutine identify_molecules(string,nmolec,molec_names)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c     I/O
      character*(Nchar_mx) string
      integer nmolec
      character*10 molec_names(1:Nmol_max)
c     temp
      character*(Nchar_mx) input_file
      integer i,j,icb,imol,iht
      logical cb_found,is_ht_mol
      character*(Nchar_mx) molname
c     HITRAN data
      integer HNmol
      integer Hmol_index(1:Nmol_max)
      integer Hmol_niso(1:Nmol_max)
      character*10 Hmol_name(1:Nmol_max)
      integer Hiso_index(1:Nmol_max,1:Niso_max)
      double precision Hiso_abundance(1:Nmol_max,1:Niso_max)
      double precision Hiso_Qref(1:Nmol_max,1:Niso_max)
      integer Hiso_g(1:Nmol_max,1:Niso_max)
      double precision Hiso_mass(1:Nmol_max,1:Niso_max)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine identify_molecules'

c     get molecules names and indexes from HITRAN
      input_file='./molparam_hitran2008.txt'
      call read_molparam(input_file,
     &     HNmol,Hmol_index,Hmol_niso,Hmol_name,
     &     Hiso_index,Hiso_abundance,Hiso_Qref,Hiso_g,Hiso_mass)
c     Debug
c      write(*,*) 'HITRAN molecules:'
c      do imol=1,HNmol
c         write(*,*) Hmol_index(imol)
c     &        ,Hmol_name(imol)(1:strlen(Hmol_name(imol)))
c      enddo ! imol
c     Debug

      nmolec=0
      do i=1,strlen(string)
         if (string(i:i).eq.'[') then
            cb_found=.false.
            do j=i+1,strlen(string)
               if (string(j:j).eq.']') then
                  cb_found=.true.
                  icb=j
                  goto 111
               endif
            enddo ! j
 111        continue
            if (.not.cb_found) then
               call error(label)
               write(*,*) 'the "]" character could not be found in:'
               write(*,*) string(i+1:strlen(string))
               stop
            endif
c     A molecule name has been located in the string
            molname=string(i+1:icb-1)
c     check if this molecule belongs to the list of HITRAN molecules
            is_ht_mol=.false.
            do imol=1,HNmol
               if (molname(1:strlen(molname)).eq.
     &              Hmol_name(imol)(1:strlen(Hmol_name(imol)))) then
                  is_ht_mol=.true.
                  iht=imol
                  goto 112
               endif
            enddo ! imol
 112        continue
            if (is_ht_mol) then
               nmolec=nmolec+1
               molec_names(nmolec)=Hmol_name(iht)
     &              (1:strlen(Hmol_name(iht)))
            else
c     Special case of HDO
               if (molname(1:strlen(molname)).eq.'HDO') then
                  nmolec=nmolec+1
                  molec_names(nmolec)='HDO'
               else
                  call error(label)
                  write(*,*) 'Molecule is not among the HITRAN list:'
                  write(*,*) molname(1:strlen(molname))
                  stop
               endif ! HDO
            endif ! is_ht_mol
         endif ! '[' has been identified
      enddo ! i

      return
      end