pro read_combo,expt,dtype,spec,year,smon,emon,mdata,lat,lon,plevs,dates,thlevs=thlevs,$
	tplevs=tplevs,dir=dir,fname=fname,reduce=reduce,trpp=trpp,nohdr=nohdr,grid=grid

;	New capability to read in metwater (H2O from amonthly files). Set spec to 'water'

;	Allow different resolutions to be read in: 4x5, 2x2.5, and 1x1.25. Default is
;	2x2.5 but res=1 is for 1x1.25 and res=4 is 4x5

;	New Oct 2007: Option to read in tropopause pressure. Currently it's available
;	 	only in idaily files.
;	New Sept 2007:  tplevs is full 4D true pressure array.
;		keyword /reduce will the frequency of idaily output to keep array sizes
;	 	down. Set reduce to the number of days per month desired. For example
;		reduce=6 retrieves 6 days per month (~5 day interval), reduce=8 
;		retrieves with ~3 day intervals (mod(28/reduce)). Dates are output.

;	February, 2007. New reader for combo output. Combo 124-species mechanism.
;		'PV' and 'temp' are also available in the daily files.
;
;	Possible datatypes (dtype) to read in are:
;	   .amonthly.nc -> monthly avg concentrations, formerly the .const.nc files
;	   .overpass1.nc -> avg concentrations 10-11am, daily
; 	   .overpass2.nc -> avg concentration 1-2pm, daily
;   	   .idaily.nc -> daily instantaneous conc at 12Z (freq1), formerly .daily.nc
;	   .hourlyoz.nc -> hourly surface ozone (freq2)
;	   .column.nc -> instantaneous daily ozone column at 12Z (freq3)
;	NOTE: HOURLY NOT YET WORKING!!!
;	
;       INPUTS: spec is a string for species name.
;               year is a 4-char string, e.g. '1995' or '2005'
;		smon is starting month, a 3-char string
;		emon is the ending month, a 3-char string	
;               expt is the experiment name (appears in file name)
;               dtype is datatype, listed above. At a minimum, type must cotain a substring
;		   of filetype that uniquely identifies it.

;	Filenames to be built are modelname_exptname_YYYY_MMM.datatype.nc

;       OPTIONAL INPUT: Set thlevs to desired theta levels to return values on theta.
;	    Define DIR if you want to specify the directory where the data is stored (if you 
;	    are not working on the GSFC code916 cluster).
;	    Set fname equal to the filename, but then you can only do one month at a time.
;	tplevs=1 if you want true, 3D pressure levels returned
;	reduce=n, where n is the number of days per month you want (if choosing a datatype
;	     that has daily output).

;	*****************************************************************
;		SET DEFAULT DIRECTORY to GMI Public Directory
;	*****************************************************************

defdir='/science/gmi/data/gmic/'

;	************** Use expt & year to determine subdirectory of files ********

exptdir=['Aura2', 'Aura2_nlp', 'Aura2_F12_nlp', 'Aura2_F24_nlp', 'G4GCM_nlp','Aura4','G5aura','ODP2000',$
    'MERRA','MERRA_noclouds','Aura4Synoz','MERRA_no_het','MERRA_gcr','MERRA_jNOp6','HindcastSpinup',$
    'Hindcast','Hindcast_noGCRs','HindcastFF','Hindcast3','MERRA-igac','Hindcast3Igac','Hindcast3Igac2',$
    'MERRA_no_hetarc','MERRA_no_hetant','MERRA_BrTest','HindcastFFIgac2','Hindcast-CCMI',$
    'HindcastFFIgac2-HighRes','HindcastFFIgac2+5Br-HighRes','HindcastFFIgac2-MO2+ClO-HighRes',$
    'HindcastFFIgac2+5Br-MO2+ClO-HighRes','HindcastFFIgac2+5Br-MO2+ClO_newrate-HighRes',$
    'HindcastFFIgac2_1JNO+5Br-MO2+ClO_newrate-HighRes','HindcastBr','HindcastFFIgac2-Merra2-HighRes',$
    'MERRA_fixODS','HindcastMR2','HindcastMR2-CCMI']

nr=n_elements(exptdir)
preslev=[42,42,42,42,42,42, fltarr(nr-6)+72]

  testdir=strlowcase(exptdir)
  lowexpt=strlowcase(expt)
  e=where(testdir eq lowexpt)
;	aura2_nlp is synonymous with aura3
  if lowexpt eq 'aura3' then e=1

if e lt 0 and keyword_set(dir) eq 0 then begin
   print, 'Possible Experiment (`expt`) Names Are:'
   print, '  `Aura2` - the original Aura run, Mar 2007'
   print, '  `Aura2_nlp or Aura3` - 2nd Aura run with New Lightning, Aug 2007'
   print, '  `G4GCM_nlp` - 2nd FVGCM (94-98) run with New Lightning, Sep 2007'
   print, '  `Aura2_F12_nlp` - Aura period, Forecast Met fields (12-24Hr Fcst)'
   print, '  `Aura2_F24_nlp` - Aura period, Forecast Met fields (12-36Hr Fcst)'
   print, '  `Aura4 - Aura period, Feb 2004-Dec 2007'
   print, '  `ODP2000 - GEOS5-GCM ODP 10-year expts'
   print, '  `G5aura - Has synoz bug, Aura period, 2005-2008.'
   print, '  `MERRA - Aura period so far.'
   print, '  `MERRA_no_het, MERRA_no_hetarc, and MERRA_no_hetant`'
   print, '  `MERRA_BrTest` '
   print, '  `MERRA_gcr - Jan 2004-2007? '
   print, '  `MERRA_jNOp6 - Jan 2004-2007? '
   print, '  `HindcastSpinup`'
   print, '  `Hindcast, HindcastFF, and HindcastFFIgac2` '
   print, '  `Hindcast_noGCRs`'
   print, '  `Hindcast3, Hindcast3Igac, and Hindcast3Igac2`'
   print, '  `Hindcast-CCMI`'
   print, '  `MERRA-igac, MERRA_fixODS`'
   print, '  `HindcastFFIgac2-HighRes, HindcastFFIgac2+5Br-HighRes, and HindcastFFIgac2-MO2+ClO-HighRes`'
   print, '  `HindcastFFIgac2+5Br-MO2+ClO_newrate-HighRes, HindcastFFIgac2_1JNO+5Br-MO2+ClO_newrate-HighRes`'
   print, '  `HindcastBr, HindcastFFIgac2-Merra2-HighRes, HindcastMR2`, HindcastMR2-CCMI'
   goto, theend
endif

;       This sets the number of pressure levels for each of these experiments
   np=preslev(e)

;	Modify default directory to include experiment and year subdirectories

if keyword_set(dir) eq 0 then begin
   dir=defdir+exptdir(e)+'/'+year
   if keyword_set(grid) eq 1 then if grid eq 1 and e eq 8 then dir=defdir+'MERRA_1x1.25/'+year
   if keyword_set(grid) eq 1 then if grid eq 1 and e eq 11 then dir=defdir+'MERRA_no_het_1x1.25/'+year
   if keyword_set(grid) eq 1 then if grid eq 1 and e eq 22 then dir=defdir+'MERRA_no_hetarc_1x1.25/'+year
   if keyword_set(grid) eq 1 then if grid eq 1 and e eq 23 then dir=defdir+'MERRA_no_hetant_1x1.25/'+year
   if keyword_set(grid) eq 1 then if grid eq 1 and e eq 24 then dir=defdir+'MERRA_BrTest_1x1.25/'+year
endif


;       **************************************************************
;       	Set Month number and number of days.
;	**************************************************************

mlab=['jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec']
ndm=[31,28,31,30,31,30,31,31,30,31,30,31]

;	Get 29 days for February in leap years
iyear=fix(year)
if iyear mod 4 eq 0 then ndm=[31,29,31,30,31,30,31,31,30,31,30,31]

smon=strlowcase(smon) & emon=strlowcase(emon)
mi=where(mlab eq smon) & mf=where(mlab eq emon)

;	******** Reject if start/end months are out of order ************

if mf(0) lt mi(0) then begin
   print, 'End month comes before Start month.'
   goto, theend
endif

;	********* Reject values of reduce lt 1 or gt 10 *********************

if keyword_set(reduce) then begin
   if reduce lt 1 or reduce gt 10 then begin
      print, 'Keyword `reduce` may have integer values between 1 and 10, and it'
      print, '  can only be used with idaily or overpass datatypes.'
      goto, theend
   endif
endif


;	***********************************************************************
;       	FIGURE OUT THE DATATYPE TO CREATE THE FILENAME
;	***********************************************************************


fdtype=['amonthly','idaily','overpass1','overpass2','overpass3','overpass4','hourlyoz','column','rst']

;		CASE where filename is input through optional keyword

if keyword_set(fname) eq 1 then begin
   nm=1
   cfile=strarr(nm)
   cfile(0)=dir+'/'+fname
   ndays=28
   subtype=['mon','dai','pass1','pass2','pass3','pass4','hour','col','rst']
   ss=n_elements(subtype)
   res=intarr(ss)
   for i=0,ss-1 do res(i)=strpos(fname,subtype(i))
   x=where(res ge 0)
   ddtype=fdtype(x(0))
   ndays=ndm(mi(0))
   if ddtype eq 'amonthly' or ddtype eq 'rst' then ndays=1

endif

;	******* CASE where filename is created from inputs to this code  *************

if keyword_set(fname) eq 0 then begin
  a=where(strmatch(fdtype,'*'+dtype+'*'))

;	*********** Reject if wrong data type used ************

  if a(0) lt 0 then begin
     print, 'Acceptable input for dtype is any unique substring of the following data types:'
     print, fdtype
     goto, theend
  endif

  ddtype=fdtype(a)
  if ddtype eq fdtype(6) then begin
     print, 'Hourly ozone reader not available.'
     goto, theend
  endif
;       

;	******** Create actual filenames for each experiment **********

   if e eq 0 then fp='/gmic_aura2_'+year+'_'
   if e eq 1 then fp='/gmic_aura2ReRunHO2Light_'+year+'_'
   if e eq 2 then fp='/gmic_aura2for12h_'+year+'_'
   if e eq 3 then fp='/gmic_aura2for24h_'+year+'_'
   if e eq 4 then fp='/gmic_fvgcm_'+year+'_'
   if e eq 5 then fp='/gmic_aura4_'+year+'_'
   if e eq 6 then fp='/gmic_G5aura_'+year+'_'
;  if e eq 7 then fp='/gmic_G5aura_'+year+'_'
   if e eq 8 then fp='/gmic_MERRA_'+year+'_'
   if e eq 9 then fp='/gmic_MERRA_noclouds_'+year+'_'
   if e eq 10 then fp='/gmic_Aura4Synoz_'+year+'_'
   if e eq 11 then fp='/gmic_MERRA_no_het_'+year+'_'
   if e eq 12 then fp='/gmic_MERRA_gcr_'+year+'_'
   if e eq 13 then fp='/gmc_MERRA_jNOp6_'+year+'_'
   if e eq 14 then fp='/gmic_HindcastSpinup_'+year+'_'
   if e eq 15 then fp='/gmic_Hindcast_'+year+'_'
   if e eq 16 then fp='/gmic_HindcastSpinupNoGCR_'+year+'_'
   if e eq 17 then fp='/gmic_HindcastFF_'+year+'_'
   if e eq 18 then fp='/gmic_Hindcast3_'+year+'_'
   if e eq 19 then fp='/gmic_MERRA-igac_'+year+'_'
   if e eq 20 then fp='/gmic_Hindcast3Igac_'+year+'_'
   if e eq 21 then fp='/gmic_Hindcast3Igac2_'+year+'_'
   if e eq 22 then fp='/gmic_MERRA1x1_no_hetarc_'+year+'_'
   if e eq 23 then fp='/gmic_MERRA1x1_no_hetant_'+year+'_'
   if e eq 24 then fp='/gmic_MERRA1x1_BrTest_'+year+'_'
   if e eq 25 then fp='/gmic_HindcastFFIgac2_'+year+'_'
   if e eq 26 then fp='/gmic_Hindcast-CCMI_'+year+'_'
   if e eq 27 then fp='/gmic_HindcastFFIgac2-HighRes_'+year+'_'
   if e eq 28 then fp='/gmic_HindcastFFIgac2+5Br-HighRes_'+year+'_'
   if e eq 29 then fp='/gmic_HindcastFFIgac2-MO2+ClO-HighRes_'+year+'_'
   if e eq 30 then fp='/gmic_HindcastFFIgac2+5Br-MO2+ClO-HighRes_'+year+'_'
   if e eq 31 then fp='/gmic_HindcastFFIgac2+5Br-MO2+ClO_newrate-HighRes_'+year+'_'
   if e eq 32 then fp='/gmic_HindcastFFIgac2_1JNO+5Br-MO2+ClO_newrate-HighRes_'+year+'_'
   if e eq 33 then fp='/gmic_HindcastBr_'+year+'_'
   if e eq 34 then fp='/gmic_HindcastFFIgac2Merra2-HighRes_'+year+'_'
   if e eq 35 then fp='/gmic_MERRA_fixODS_'+year+'_'
   if e eq 36 then fp='/gmic_HindcastMR2_'+year+'_'
   if e eq 37 then fp='/gmic_HindcastMR2-CCMI_'+year+'_'


;	The original 1x1.25 run with emissions problems had the file name fp='/gmic_g5-1x1_'+year+'_'
;	The correct 1x1.25 run has the name 'MERRA_1x1'

if keyword_set(grid) eq 1 then if e eq 8  and grid eq 1 then fp='/gmic_MERRA1x1_'+year+'_'
if keyword_set(grid) eq 1 then if e eq 11 and grid eq 1 then fp='/gmic_MERRA1x1_'+year+'_'

   nm=mf(0)-mi(0)+1
   cfile=strarr(nm)
   ndays=0
   for t=mi(0),mf(0) do ndays=ndays+ndm(t)
   for m=0,nm-1 do cfile(m)=dir+fp+mlab(mi(0)+m)+'.'+ddtype+'.nc'
   if ddtype eq 'amonthly' then ndays=nm

endif


;       ***************************************************************************
;       Convert input 'spec' string to species number. Determine variable name for filetype.
;	***************************************************************************
;
specname=spec

;       New capability to read in variable 'metwater'. It can be in idaily or amonthly

if spec eq 'water' and dtype eq 'amonthly' then begin
    varname='metwater'
    ddtype='amonthly'
    goto, skipvarname
endif

if spec eq 'water' and dtype eq 'idaily' then begin
    varname='metwater_freq1'
    ddtype='idaily'
    goto, skipvarname
endif

;	Now using convert_species to replace convert_spec124_new. This works for 119 and 124 species mechanisms.

spec=strupcase(spec)
convert_species,cfile(0),ddtype,spec,spnum
if spnum eq 0 then begin
   print, 'Species is a string.'
   print, 'The species you chose, ',spec,', is either not in the combo mechanism,'
   print, '   or is not saved in this filetype (e.g., not saved in overpass1).'
   print, 'Type @get_species_combo124 to see all possible species.'
   print, 'PV and Temp are also species.'
   goto, theend
endif

varname=findvarname(ddtype,spnum,spec)

skipvarname: 

;	***********************************************************************
;			READ IN Default 2x2.5 (latxlon) MODEL DATA or change resolution
;	***********************************************************************

if keyword_set(thlevs) then nk=n_elements(thlevs) else nk=np

ni=144 & nj=91
if keyword_set(grid) eq 1 then begin
   if grid eq 1 then ni=288
   if grid eq 1 then nj=181
   if grid eq 4 then ni=72
   if grid eq 4 then nj=46
endif

mdata=fltarr(ni,nj,nk,ndays) 
xplevs=fltarr(ni,nj,np,ndays)
xtrpp=fltarr(ni,nj,ndays)

if keyword_set(reduce) then begin
   mdata=fltarr(ni,nj,nk,nm*reduce)
   xplevs=fltarr(ni,nj,np,nm*reduce)
   xtrpp=fltarr(ni,nj,nm*reduce)
endif

if ddtype eq 'idaily' or ddtype eq 'overpass1' or ddtype eq 'overpass2' or $
  ddtype eq 'overpass3' or ddtype eq 'overpass4' then begin
  lastn=0
  for m=0,nm-1 do begin 
     confile=cfile(m) 
     nrec=ndm(mi(0)+m)
     ind=indgen(ndm(mi(0)+m))+lastn
     print, 'Reading ',specname,' from ',confile
     read_const,confile,spnum,nrec,lat,lon,plevs,data,thlevs=thlevs,type=varname,tplevs=tplevs,$
        ddtype=ddtype,reduce=reduce,trpp=trpp,nohdr=nohdr
;	if read_const is returning a constituent, it will be scaled to ppb. PV and temp unaffected

     if keyword_set(reduce) then ind=indgen(reduce)+lastn
     mdata(*,*,*,ind)=data 
     if keyword_set(tplevs) and keyword_set(thlevs) eq 0 then xplevs(*,*,*,ind)=tplevs
     if keyword_set(trpp) then xtrpp(*,*,ind)=trpp
     lastn=max(ind)+1
  endfor
endif

if ddtype eq 'amonthly' or ddtype eq 'rst' then begin
  for m=0,nm-1 do begin
     confile=cfile(m)
     print, 'Reading ',specname,' from ',confile
     read_const,confile,spnum,1,lat,lon,plevs,data,thlevs=thlevs,type=varname,tplevs=tplevs,ddtype=ddtype,nohdr=nohdr
     mdata(*,*,*,m)=data
     if keyword_set(tplevs) then xplevs(*,*,*,m)=tplevs
  endfor
endif

if ddtype eq 'column' then begin
   mdata=fltarr(ni,nj,ndays)
   if keyword_set(reduce) then mdata=fltarr(ni,nj,nm*reduce)
   lastn=0
   for m=0,nm-1 do begin
     confile=cfile(m)
     nrec=ndm(mi(0)+m)
     ind=indgen(ndm(mi(0)+m))+lastn
     print, 'Reading ',varname,' ',specname,' from ',confile
     read_const,confile,spnum,nrec,lat,lon,plevs,data,thlevs=thlevs,type=varname,ddtype=ddtype,reduce=reduce
;
     if keyword_set(reduce) then ind=indgen(reduce)+lastn
     mdata(*,*,ind)=data
     lastn=max(ind)+1
  endfor
endif


if ddtype ne 'column' then print, 'Constituents are scaled to ppb.'

if keyword_set(tplevs) and keyword_set(thlevs) eq 0 then tplevs=xplevs
trpp=xtrpp

;	*************************************************************
;			Create array of Dates 
;	*************************************************************

mm=['01','02','03','04','05','06','07','08','09','10','11','12']
dd=[mm,strtrim(string(indgen(19)+13),2)]

dates=strarr(nm)
if ddtype eq 'amonthly' or ddtype eq 'rst' then for i=0,nm-1 do dates(i)=year+mm(mi(0)+i)
if ddtype ne 'amonthly' and ddtype ne 'rst' and keyword_set(reduce) eq 0 then begin
   dates=strarr(ndays)
   ind=0
   for i=mi(0),mf(0) do begin
      for n=0,ndm(i)-1 do dates(ind+n)=year+mm(i)+dd(n)
      ind=ind+ndm(i)
   endfor
endif

if ddtype ne 'amonthly' and ddtype ne 'rst' and keyword_set(reduce) then begin
   dates=strarr(ndays)

;	Figure out which days per month were kept. This is directly from read_const

   if reduce ge 5 then begin
          int=fix(27/(reduce-1))
          rx=indgen(reduce)*int
       endif

   if reduce lt 5 then begin
          if reduce eq 1 then rx=[0]
          if reduce eq 2 then rx=[0,14]
          if reduce eq 3 then rx=[0,10,20]
          if reduce eq 4 then rx=[0,7,14,21]
       endif

   nrx=n_elements(rx)
   ind=0
   dates=strarr(nrx*(mf(0)-mi(0)+1))
   for i=mi(0),mf(0) do begin
      for n=0,nrx-1 do dates(ind+n)=year+mm(i)+dd(rx(n))
      ind=ind+nrx
   endfor
endif

theend:
end