pro read_column_multi,expt,spec,syear,eyear,mdata,mmdata,lat,lon,ylab,grid=grid,reduce=reduce,zm=zm

;	zm=1 returns zonal mean.
;	Add reduced number of days per month. reduce must be an integer between 1-10.
;	Read in multiple years of column output (lon,lat,days,years)
;	Calls read_combo ny times, where ny is number of years. Option grid for reading
;	in 1x1.25 runs.

leaplist=indgen(15)*4+1976

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

;	Create an array matching the years that identifies leap years.

for y=0,ny-1 do if iyear(y) mod 4 eq 0 then lymask(y)=1

days=[indgen(59),indgen(306)+60]

;	Create output array
ni=144 & nj=91 
if keyword_set(grid) then ni=288
if keyword_set(grid) then nj=181

if keyword_set(zm) eq 1 then mdata=fltarr(nj,365,ny) else mdata=fltarr(ni,nj,365,ny)
if keyword_set(reduce) then mdata=fltarr(ni,nj,reduce*12,ny)
if keyword_set(reduce) and keyword_set(zm) eq 1 then mdata=fltarr(nj,reduce*12,ny)

for y=0,ny-1 do begin 
  read_combo,expt,'column',spec,ylab(y),'jan','dec',data,lat,lon,hplevs,dates,grid=grid,reduce=reduce
  if keyword_set(zm) eq 1 then begin
     if lymask(y) gt 0 and keyword_set(reduce) eq 0 then mdata(*,*,y)=avg(data(*,*,days),0) $
          else mdata(*,*,y)=avg(data,0)
  endif
  if keyword_set(zm) eq 0 then begin
     if lymask(y) gt 0 and keyword_set(reduce) eq 0 then mdata(*,*,*,y)=data(*,*,days) else mdata(*,*,*,y)=data
  endif

endfor

;	Create monthly averaged output too when reduce is not used

ndm=[31,28,31,30,31,30,31,31,30,31,30,31]
di=intarr(12) & df=di
for i=1,11 do di(i)=total(ndm(0:i-1))
for i=0,11 do df(i)=di(i)+ndm(i)-1

if zm eq 0 then mmdata=fltarr(ni,nj,12,ny) else mmdata=fltarr(nj,12,ny)
if keyword_set(reduce) eq 0 and keyword_set(zm) eq 0 then $
   for y=0,ny-1 do for m=0,11 do mmdata(*,*,m,y)=avg(mdata(*,*,di(m):df(m),y),2)
if keyword_set(reduce) eq 0 and keyword_set(zm) eq 1 then $
   for y=0,ny-1 do for m=0,11 do mmdata(*,m,y)=avg(mdata(*,di(m):df(m),y),1)


end