pro usage
  print,format="(/,72('-'))"
  print,format="('readfile.pro contains functions and procedures to read a TGCM history file.')"
  print,format="(/,'Following are examples of usage:',/)"
  print,'  h = readfile(/help)     ; get help on how to use function readfile'
  print,'  h = getvar(/help)       ; get help on how to use function getvar'
  print,'  showvar,/help           ; get help on how to call procedure showvar'
  print,' '
  print,'  file = readfile()       ; a pickfile dialogue box will pop up'
  print,'  help,file,/struct       ; get help on file structure (after readfile)'
  print,'  print,file.f4dnames     ; print names of 4d fields on file
  print,'  print,file.varnames     ; print names of all variables on the file
  print,'  h = getvar(/help)       ; get help using getvar
  print,'  TN = getvar(file.ncid,''TN'',tn_info)   ; read variable TN (and get info structure)'
  print,'  help,tn_info,/struct    ; get help on returned structure tn_info'
  print,'  print,min(TN),max(TN)   ; print min,max of TN'
  print,' ' 
  print,'  file = readfile(''/toshi/ftp/pub/tgcm/data'') ; a pickfile dialogue box will pop up'
  print,'  showvar,/help           ; get help using showvar'
  print,'  showvar,file.ncid,''TN''  ; print info about field TN'
  print,format="(72('-'),/)"
end
;-----------------------------------------------------------------------
function openfile,file
;
; Open netcdf file. If successful, return cdfid,
; otherwise return -1.
;
if (strlen(strcompress(file,/remove_all)) le 0) then begin
  print,'>>> WARNING openfile: empty file name'
  return,-1
endif
on_ioerror,badfile
cdfid = ncdf_open(file)
print,' '
print,'Opened file ',file
return,cdfid
badfile:
print,' '
print,'>>> openfile: error opening file ',file
print,' '
return,-1
end
;-----------------------------------------------------------------------
function abort,path
  file = { path:' ' }
  file.path = path
  return,file
end
;-----------------------------------------------------------------------
function getdim,ncid,name
id = ncdf_dimid(ncid,name)
if (id gt -1) then begin
  ncdf_diminq,ncid,id,name,len
  return,len
endif else begin
  print,'>>> WARNING: could not find "',name,'" dimension.'
  len = 0
  return,len
endelse
end
;-----------------------------------------------------------------------
function getvar,ncid,name,varinfo,help=help
if (keyword_set(help)) then begin
  print,format="(/,72('-'))"
  print,format="('getvar: read a variable from a currently open file')"
  print,format="(/,'Usage (called after a call to readfile)')"
  print,format="('  var = getvar(ncid,name,varinfo,[/help])')"
  print,format="(/,'Where:')"
  print,format="('  ncid is the current file id')"
  print,format="('  name is the name of the variable to read')"
  print,format="('  varinfo is a returned structure containing useful information about  the variable')"
  print,format="('  help is an optional keyword, to get this useage message.')"
  print,format="(/,'Example 1: Get this useage message')"
  print,format="('  h = getvar(/help)')"
  print,format="(/,'Example 2 (after calling readfile to get file structure): Read field TN')"
  print,'  TN = getvar(file.ncid,''TN'',vinfo)'
  print,' ' 
  print,'After the above call, execute ''help,vinfo,/struct'' to view the contents'
  print,'  of the returned structure.'
  print,format="(72('-'),/)"
  return,abort('')
endif
id = ncdf_varid(ncid,name)
if (id gt -1) then begin
; print,'Reading variable ',name,'...'
  ncdf_varget,ncid,id,f
  varinfo = ncdf_varinq(ncid,id)
  return,f
endif else begin
  print,'>>> WARNING: could not find variable with name "',name,'".'
  return,0
endelse
end
;-----------------------------------------------------------------------
pro showvar,ncid,name,help=help
;
; Usage statement:
;
if (keyword_set(help)) then begin
  print,format="(/,72('-'))"
  print,format="('showvar: print useful information about a variable on the current file')"
  print,format="(/,'Usage (usually called after a call to readfile):')
  print,format="('  showvar,ncid,varname')"
  print,format="(/,'Where:')"
  print,format="('  ncid is the current file id')"
  print,format="('  varname is the name of the desired variable')" 
  print,format="(/,'Example 1 (after reading file structure): Show field TN')"
  print,'  showvar,file.ncid,''TN'''
  print,format="(/,'Example 2: Show 1st 4d field')"
  print,format="('  showvar,file.ncid,file.f4dnames[0]')"
  print,format="(72('-'),/)"
  return
endif
id = ncdf_varid(ncid,name)
if (id gt -1) then begin
  varinfo = ncdf_varinq(ncid,id)
; help,/str,varinfo
  print,format="(/,'Field ',a,':')",name
  print,format="('  Datatype: ',a)",varinfo.datatype
  if (varinfo.ndims gt 0) then begin
    dimnames = strarr(varinfo.ndims)
    dimlens = strarr(varinfo.ndims)
    print,format="('  There are ',i3,' dimensions: ')",varinfo.ndims
    for i=0,varinfo.ndims-1 do begin
      ncdf_diminq,ncid,varinfo.dim[i],dimname,len
      dimnames[i] = dimname 
      dimlens[i] = len 
      print,format="('    ',a,' = ',i4)",dimnames[i],dimlens[i]
    endfor
  endif
  if (varinfo.natts gt 0) then begin
    print,format="('  There are ',i3,' attributes: ')",varinfo.natts
    for i=0,varinfo.natts-1 do begin
      attname = ncdf_attname(ncid,id,i)
      att = ncdf_attinq(ncid,id,attname)
;     help,/str,att
      ncdf_attget,ncid,id,attname,value
      if (att.datatype eq 'CHAR') then $
        print,format="('    ',a,' = ',a)",attname,string(value)
    endfor
  endif
endif else begin
  print,'>>> WARNING: could not find variable with name "',name,'".'
endelse
print,' '
end
;-----------------------------------------------------------------------
function readfile,path=path,help=help
;
; Usage statement:
;
if (keyword_set(help)) then begin
  print,format="(/,72('-'))"
  print,format="('readfile: Function to read tiegcm or timegcm history file, and to')"
  print,format="('          return a data structure containing useful information.')"
  print,format="(/,'Usage:')"
  print,format="('  result = READFILE([PATH=string][,/HELP])')"
  print,format="(/,'Where:')"
  print,'  ''string'' is an optional path to a directory or file.'
  print,'  ''result'' is the returned structure describing the file.'
  print,format="(/,'Example 1: Get this usage message:')"
  print,format="('  help = readfile(/help)')"
  print,format="(/,'Example 2: Call without a path (a pickfile dialogue will pop up)')"
  print,format="('  file = readfile()')"
  print,format="(/,'Example 3: Call with a directory (a pickfile dialogue will pop up)')"
  print,'  file = readfile(path=''/home/mylogin/dir1/dir2'')'
  print,format="(/,'Example 4: Read file tiegcm.nc')"
  print,'  file = readfile(path=''tiegcm.nc'')'
  print,' ' & print,'After any of the above calls, execute ''help,file,/struct'' to view the contents of the structure.'
  print,format="(72('-'),/)"
  return,abort('')
endif
;
; Path was provided.
; If path does not exist, or is not readable, abort with error.
; If path is a directory, pass it to dialog_pickfile
;
if (keyword_set(path)) then begin
  fileinfo = file_info(path)
  if (not fileinfo.exists) then begin
    print,format="(/,'>>> readfile: file or directory does not exist: ',a,/)",path
    return,abort(path)
  endif else if (not fileinfo.read) then begin
    print,format="(/,'>>> readfile: file ',a,' is not readable.')",path
    return,abort(path)
  endif else if (fileinfo.directory) then begin
    title = 'Please choose a netCDF history file'
    filepath = dialog_pickfile(title=title,path=path,/read,/must_exist,filter='*.nc')
  endif else begin
    filepath = path 
  endelse
endif else begin ; path not provided
  ;path = '/hao/tgcm/timegcm1.2_yr04/sy4n'
  ;path = '/toshi/ftp/pub/foster/murray'
   path = '/e/foster/tiegcm1.8/tiegcm1.8-linux/FOSTER.tiegcm1.8.stest001.nc'
  ;path = ' '
  if (strlen(strcompress(path,/remove_all)) le 0) then begin
    tgcmdata = getenv('TGCMDATA')
    if (strlen(strcompress(tgcmdata,/remove_all)) gt 0) then begin
      path = tgcmdata
    endif else begin
      cd,current=path
    endelse
  endif
  title = 'Please choose a netCDF history file'
  filepath = dialog_pickfile(title=title,path=path,/read,/must_exist,filter='*.nc')
endelse
;
; Open history file:
;
ncid = openfile(filepath)
if (ncid lt 0) then return,abort(filepath)
ncdf_control,ncid,/verbose
;
fileinfo = ncdf_inquire(ncid)
;help,fileinfo,/struct
;
; Get grid dimensions (lat,lon,lev and time are required):
;
; Geographic dimensions:
nlat = getdim(ncid,'lat') & if (nlat lt 0) then return,abort(filepath)
nlon = getdim(ncid,'lon') & if (nlon lt 0) then return,abort(filepath)
nlev = getdim(ncid,'lev') & if (nlev lt 0) then return,abort(filepath)
ntime = getdim(ncid,'time') & if (ntime lt 0) then return,abort(filepath)
nilev = getdim(ncid,'ilev')
;
; Magnetic dimensions:
nmlat = getdim(ncid,'mlat')
nmlon = getdim(ncid,'mlon')
nmlev = getdim(ncid,'mlev')
nimlev = getdim(ncid,'imlev')
;
;print,format="('nlat=',i4,' nlon=',i4,' nlev=',i4,' nilev=',i4,' ntime=',i4)",$
;  nlat,nlon,nlev,nilev,ntime
;
; Get coordinate arrays:
;
; Geographic coords:
lat = getvar(ncid,'lat') & ;print,'lat=',lat
lon = getvar(ncid,'lon') & ;print,'lon=',lon
lev = getvar(ncid,'lev') & ;print,'lev=',lev
if (nilev gt 0) then ilev = getvar(ncid,'ilev') else ilev = 0
;
; Magnetic coords:
if (nmlat gt 0) then mlat = getvar(ncid,'mlat') else mlat = 0
if (nmlon gt 0) then mlon = getvar(ncid,'mlon') else mlon = 0
if (nmlev gt 0) then mlev = getvar(ncid,'mlev') else mlev = 0
if (nimlev gt 0) then imlev = getvar(ncid,'imlev') else imlev = 0
time = getvar(ncid,'time') & ;print,'time=',time
;
; Get model times:
;
mtime = getvar(ncid,'mtime') & ;print,'mtime=',mtime
;
; Get list of 4-d vars:
;
n4d = 0
for i=0,fileinfo.nvars-1 do begin
  varinfo = ncdf_varinq(ncid,i)
  if (varinfo.ndims eq 4) then n4d = n4d+1
endfor
f4dnames = strarr(n4d)
varnames = strarr(fileinfo.nvars)
n4d = 0
for i=0,fileinfo.nvars-1 do begin
  varinfo = ncdf_varinq(ncid,i)
; help,/str,varinfo
  varnames[i] = varinfo.name
  if (varinfo.ndims eq 4) then begin
    n4d = n4d+1
    f4dnames(n4d-1) = varinfo.name
;   f = getvar(ncid,varinfo.name,vinfo)  ; this read the 4d var
  endif
endfor
;
; Return file structure:
;
file = {        	$
  path:filepath,	$
  ncid:ncid,    	$
  ntime:ntime,  	$
  mtime:mtime,  	$
  lat:lat,      	$
  lon:lon,      	$
  lev:lev,      	$
  ilev:ilev,    	$
  mlat:mlat,            $
  mlon:mlon,            $
  mlev:mlev,            $
  imlev:imlev,          $
  f4dnames:f4dnames,	$
  varnames:varnames 	$
}
;
; Report to stdout:
;
;help,file,/struct
print,format="(/,72('-'))"
print,format="('File ',a)",file.path
print,format="('  ntime=',i4)",file.ntime
;print,format="('  mtime=',(3i4))" & print,file.mtime
print,format="('  First, last mtime = ',3i4,',',3i4)",$
  file.mtime[*,0],file.mtime[*,file.ntime-1]
print,format="('  nlon =',i4,' nlat =',i4,' nlev =',i4,' nilev =',i4)",$
  nlon,nlat,nlev,nilev
print,format="('  nmlon=',i4,' nmlat=',i4,' nmlev=',i4,' nimlev=',i4)",$
  nmlon,nmlat,nmlev,nimlev

print,format="(/,'List of all variables:',/)" & print,file.varnames
print,format="(/,'List of 4d fields:',/)" & print,file.f4dnames

print,format="(72('-'),/)"
return,file     ; return file structure
end             ; function readfile

