subroutine su_gases

  use gases_h

  implicit none

  integer igas, ierr, count

  !==================================================================
  !     
  !     Purpose
  !     -------
  !     Load atmospheric composition info
  !     
  !     Authors
  !     -------
  !     R. Wordsworth (2011)
  !     Allocatable arrays by A. Spiga (2011)
  !     
  !==================================================================
  

!$OMP MASTER
  ! load gas names from file 'gases.def'
  open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
  if (ierr.eq.0) then
     write(*,*) "sugases.F90: reading file gases.def"
     read(90,*)
     read(90,*,iostat=ierr) ngasmx
     if (ierr.ne.0) then
        write(*,*) "su_gases.F90: error reading number of gases"
        write(*,*) "   (first line of gases.def) "
        call abort_physic("gases.def", "Error reading number of gases", 1)
     endif

     print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..."

     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
     do igas=1,ngasmx
        read(90,*,iostat=ierr) gnom(igas)
        if (ierr.ne.0) then
           write(*,*) 'su_gases.F90: error reading gas names in gases.def...'
           call abort_physic("gases.def", "Error reading gas names", 1)
        endif
     enddo                  !of do igas=1,ngasmx

     vgas=0
     if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))
     do igas=1,ngasmx
        read(90,*,iostat=ierr) gfrac(igas)
        if (ierr.ne.0) then
           write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
           call abort_physic("gases.def", "Error reading gas molar fractions", 1)
        endif

        ! find variable gas (if any)
        if(gfrac(igas).eq.-1.0)then 
           if(vgas.eq.0)then
              vgas=igas
           else
              print*,'You seem to be choosing two variable gases'
              print*,'Check that gases.def is correct'
              call abort_physic("gases.def", "There is at least two variable gases, but only one is allowed", 1)
           endif
        endif

     enddo                  !of do igas=1,ngasmx


     ! assign the 'igas_X' labels
     count=0
     IF (.NOT.ALLOCATED(massmol))           ALLOCATE(massmol(ngasmx))
     do igas=1,ngasmx
        if (trim(gnom(igas)).eq."H2_" .or. trim(gnom(igas)).eq."H2") then
           igas_H2=igas
           massmol(igas_H2)=2.016
           count=count+1
        elseif (trim(gnom(igas)).eq."He_" .or. trim(gnom(igas)).eq."He") then
           igas_He=igas
           massmol(igas_He)=4.002602
           count=count+1
        elseif (trim(gnom(igas)).eq."H2O") then
           igas_H2O=igas
           massmol(igas_H2O)=18.01528
           count=count+1
        elseif (trim(gnom(igas)).eq."CO2") then
           igas_CO2=igas
           massmol(igas_CO2)=44.010
           count=count+1
        elseif (trim(gnom(igas)).eq."CO_" .or. trim(gnom(igas)).eq."CO") then
           igas_CO=igas
           massmol(igas_CO)=28.010
           count=count+1
        elseif (trim(gnom(igas)).eq."N2_" .or. trim(gnom(igas)).eq."N2") then
           igas_N2=igas
           massmol(igas_N2)=28.0134
           count=count+1
        elseif (trim(gnom(igas)).eq."O2_" .or. trim(gnom(igas)).eq."O2") then
           igas_O2=igas
           massmol(igas_O2)=31.998
           count=count+1
        elseif (trim(gnom(igas)).eq."SO2") then
           igas_SO2=igas
           massmol(igas_SO2)=64.066
           count=count+1
        elseif (trim(gnom(igas)).eq."H2S") then
           igas_H2S=igas
           massmol(igas_H2S)=34.082
           count=count+1
        elseif (trim(gnom(igas)).eq."CH4") then
           igas_CH4=igas
           massmol(igas_CH4)=16.0425
           count=count+1
        elseif (trim(gnom(igas)).eq."NH3") then
           igas_NH3=igas
           massmol(igas_NH3)=17.031
           count=count+1
        elseif (trim(gnom(igas)).eq."C2H6") then
           igas_C2H6=igas
           massmol(igas_C2H6)=30.069
           count=count+1
        elseif (trim(gnom(igas)).eq."C2H2") then
           igas_C2H2=igas
           massmol(igas_C2H2)=26.038
           count=count+1
        elseif (trim(gnom(igas)).eq."Ar") then
           igas_Ar=igas
           massmol(igas_Ar)=39.948
           count=count+1
        !! GG MODIF 11 Jan 2019
       elseif (trim(gnom(igas)).eq."OCS") then
           igas_OCS=igas
           massmol(igas_OCS)=60.075
           count=count+1
       elseif (trim(gnom(igas)).eq."HCl") then
           igas_HCl=igas
           massmol(igas_HCl)=36.461
           count=count+1
       elseif (trim(gnom(igas)).eq."HF") then
           igas_HF=igas
           massmol(igas_HF)=20.00634
           count=count+1
        endif
     enddo

     if(count.ne.ngasmx)then
        print*,'Mismatch between ngas and number of recognised gases in sugas_corrk.F90.'
        print*,'Either we haven`t managed to assign all the gases, or there are duplicates.'
        print*,'Please try again.'
     endif

  else
     call abort_physic("su_gases", " Unable to find gases.def", 1)
  endif
  close(90)
!$OMP END MASTER
!$OMP BARRIER

end subroutine su_gases