pro read_combo_multi,expt,spec,syear,eyear,mdata,lat,lon,hplevs,ylab,thlevs=thlevs,zm=zm,tplevs=tplevs,ht=ht,grid=grid


;	Now can read high res (grid=1) output, but only for zonal mean!
;	New optional keyword to return limited set of pressure levels. ht is an integer array
;	up to length 72 (current number of pressure levels). E.g., ht=indgen(12)+19 returns levels
;	19-30 (approx 618-192 mb).

;	Default is to return 3D fields (lon,lat,lev), but set zm=1 for zonal mean output.
;	Calls read_combo ny times, where ny is number of years. For now only do monthly means.

ny=fix(eyear)-fix(syear)+1
iyear=indgen(ny)+fix(syear)
ylab=strtrim(string(iyear),2)

;	read in first year, then get array dimensions

print, spec
if keyword_set(tplevs) then truepres=1 else truepres=0

;	For some reason read_combo is returning tplevs even when the keyword isn't set. ???
read_combo,expt,'amonthly',spec,syear,'jan','jan',data,lat,lon,hplevs,dates,thlevs=thlevs,tplevs=tplevs,grid=grid
s=size(data)
ni=s(1) & nj=s(2) & nk=s(3)
if keyword_set(ht) then nk=n_elements(ht)

if truepres eq 1 then allplevs=fltarr(ni,nj,nk,12,ny)
if truepres eq 1 and keyword_set(zm) then allplevs=fltarr(nj,nk,12,ny)

mdata=fltarr(ni,nj,nk,12,ny)
if keyword_set(zm) then mdata=fltarr(nj,nk,12,ny)

;	Don't force high res output to be zonal mean!
;if keyword_set(grid) then begin
;  mdata=fltarr(nj,nk,12,ny)
;  zm=1
;endif

if truepres eq 0 then tplevs=0

for y=0,ny-1 do begin 
  read_combo,expt,'amonthly',spec,ylab(y),'jan','dec',data,lat,lon,hplevs,dates,thlevs=thlevs,tplevs=tplevs,grid=grid
  if keyword_set(ht) then begin
     if keyword_set(zm) then mdata(*,*,*,y)=avg(data(*,*,ht,*),0) else mdata(*,*,*,*,y)=data(*,*,ht,*)
     if truepres eq 1 and keyword_set(zm) eq 0 then allplevs(*,*,*,*,y)=tplevs(*,*,ht,*) 
     if truepres eq 1 and keyword_set(zm) eq 1 then allplevs(*,*,*,*,y)=avg(tplevs(*,*,ht,*),0)
  endif
;
  if keyword_set(ht) eq 0 then begin
     if keyword_set(zm) then mdata(*,*,*,y)=avg(data,0) else mdata(*,*,*,*,y)=data
     if truepres eq 1 and keyword_set(zm) eq 0 then allplevs(*,*,*,*,0)=tplevs
     if truepres eq 1 and keyword_set(zm) eq 1 then allplevs(*,*,*,*,0)=avg(tplevs,0)
  endif

endfor

if truepres eq 1 then tplevs=allplevs

end