pro rdmet,vmet,dir,time_res,time,qnt, $
    met

; Read geos met output. 1 file per day or 1 file per month.
; INPUT (all characters):
;   vmet: geos4 or geos5
;   dir:  met directory (full path, ending with '/')
;   time_res: 3HINST,1DAVG,1MAVG
;   time: yyyymmdd (for 1D or 3H) or yyyymm (for 1M) time resolution files
;   qnt: varialbe names saved in the met files (e.g. area, airvol), plus
;        additional variables which require calculation (e.g. delz, alt, press)
; OUTPUT:
;   met: requested met fields 

fn=dir+vmet+'_met_'+time_res+'_'+time+'.nc'
cdfid=ncdf_open(fn)
if qnt eq 'press' or qnt eq 'ipress' or qnt eq 'delz' or qnt eq 'alt' $
   or qnt eq 'ialt' then begin
 nlev=72
 varname='lon'
 id=ncdf_varid(cdfid,varname)
 ncdf_varget,cdfid,id,var
 nlon=n_elements(var)
 varname='lat'
 id=ncdf_varid(cdfid,varname)
 ncdf_varget,cdfid,id,var
 nlat=n_elements(var)
 varname='area'
 id=ncdf_varid(cdfid,varname)
 ncdf_varget,cdfid,id,var
 area=var
 if qnt eq 'press' or qnt eq 'ipress' then begin
  p0=1013.25
  varname='ps'
  id=ncdf_varid(cdfid,varname)
  ncdf_varget,cdfid,id,var
  ps=var
  if qnt eq 'press' then begin
   met=fltarr(nlon,nlat,nlev)
   varname='hyam'
   id=ncdf_varid(cdfid,varname)
   ncdf_varget,cdfid,id,var
   a=var
   varname='hybm'
   id=ncdf_varid(cdfid,varname)
   ncdf_varget,cdfid,id,var
   b=var
  endif
  if qnt eq 'ipress' then begin
   met=fltarr(nlon,nlat,nlev+1)
   varname='hyai'
   id=ncdf_varid(cdfid,varname)
   ncdf_varget,cdfid,id,var
   a=var
   varname='hybi'
   id=ncdf_varid(cdfid,varname)
   ncdf_varget,cdfid,id,var
   b=var
  endif
  for i=0,nlon-1 do begin
  for j=0,nlat-1 do begin
  for l=0,nlev-1 do begin
   met(i,j,l)=a(l)*p0+b(l)*ps(i,j)
  endfor
  endfor
  endfor
 endif
 if qnt eq 'delz' or qnt eq 'alt' or qnt eq 'ialt' then begin
  varname='airvol'
  id=ncdf_varid(cdfid,varname)
  ncdf_varget,cdfid,id,var
  airvol=var
  varname='area'
  id=ncdf_varid(cdfid,varname)
  ncdf_varget,cdfid,id,var
  area=var
  delz=airvol
  for l=0,nlev-1 do delz(*,*,l)=airvol(*,*,l)/area(*,*)
  if qnt eq 'delz' then met=delz
  if qnt eq 'alt' or qnt eq 'ialt' then begin
   alt=fltarr(nlon,nlat,nlev)
   ialt=fltarr(nlon,nlat,nlev+1)
   for i=0,nlon-1 do begin
   for j=0,nlat-1 do begin
    ze=fltarr(nlev+1)
    zc=fltarr(nlev)
    for l=1,nlev do ze(l)=ze(l-1)+delz(i,j,l-1)
    for l=0,nlev-1 do zc(l)=(ze(l)+ze(l+1))/2.
    alt(i,j,*)=zc(*)
    ialt(i,j,*)=ze(*)
   endfor
   endfor
   if qnt eq 'alt' then met=alt 
   if qnt eq 'ialt' then met=ialt 
  endif
 endif 
endif else begin
 varname=qnt
 id=ncdf_varid(cdfid,varname)
 ncdf_varget,cdfid,id,var
 ncdf_close,cdfid
 met=var
endelse

end