;
function getvars,info
;
; Find variables on cdfid whose dimension names include "lat", "lon",
; and "time". This should include all 3d and 4d vars (with or without
; a vertical dimension) on both geographic and magnetic grids. 
; These vars are saved in an array of field structures 
;
printfield = 0
cdfid = info.cdfid
ncdf_control,cdfid,/verbose
file_inq = ncdf_inquire(cdfid)
nvars = 0
for ivar = 0,file_inq.nvars-1 do begin
  var = ncdf_varinq(cdfid,ivar)
  if (var.ndims gt 0) then begin
    dimnames = ' '
    for idim = 0,var.ndims-1 do begin
      ncdf_diminq,cdfid,var.dim(idim),dimname,dimsize
      dimnames = dimnames + ' ' + dimname
    endfor
;
; Group files that have (time,lat,lon) in their dimension names:
;
    time_substr = strpos(dimnames,'time')
    lat_substr = strpos(dimnames,'lat')
    lon_substr = strpos(dimnames,'lon')
    if (time_substr ge 0 and lat_substr ge 0 and lon_substr ge 0) then begin
      nvars = nvars+1        
    endif
  endif
endfor ; ivar=0,nvars-1
;
; Define a structure to hold variables from the file:
;
; Structure for a variable with lon,lat,time, and possibly lev dimensions.
; Grid may be geographic (nlon,nlat, etc), or magnetic (mlon,mlat, etc)
; (3d field if nlev==0, otherwise 4d field)
;
field = {fieldstruct,           $ ; a field structure type
  name:' ',                     $ ; short name (char8)
  long_name: ' ',               $ ; long name (may be blank)
  units: ' ',                   $ ; units (may be blank)
  idvar:0,                      $ ; netcdf variable id
  nlon:0,                       $ ; longitude dimension size
  nlat:0,                       $ ; latitude dimension size
  nlev:0,                       $ ; vertical dimension size (may be 0)
  ntime:0,                      $ ; time dimension size
  grid_type:'',                 $ ; geographic or magnetic
  dimnames:strarr(4),           $ ; dimension names in order (4th may be blank)
  ncfile:'',                    $ ; file from which field is read
  difftype:'',                  $ ; difference type ('' or '%')
  lons:ptr_new(/allocate_heap), $ ; longitude coord var
  lats:ptr_new(/allocate_heap), $ ; latitude coord var
  levs:ptr_new(/allocate_heap), $ ; levels coord var
  levname:'',levunits:'',       $ ; name and units of lev coord var (zp or height)
  levshortname:'',              $ ; short name of lev coord var
  times:ptr_new(/allocate_heap),$ ; times coord var
  data:ptr_new(/allocate_heap)  $ ; the data
}
;
; Make array of field structures:
;
if (nvars gt 0) then begin
  fields = replicate({fieldstruct},nvars)
  n = -1
  for ivar = 0,file_inq.nvars-1 do begin
    var = ncdf_varinq(cdfid,ivar)
    if (var.ndims gt 0) then begin
      dimnames = ' '
      for idim = 0,var.ndims-1 do begin
        ncdf_diminq,cdfid,var.dim(idim),dimname,dimsize
        dimnames = dimnames + ' ' + dimname
      endfor
      time_substr = strpos(dimnames,'time')
      lat_substr = strpos(dimnames,'lat')
      lon_substr = strpos(dimnames,'lon')
      if (time_substr ge 0 and lat_substr ge 0 and lon_substr ge 0) then begin
        n = n+1
;
; Set name and netcdf var id:
;
        fields[n].name = var.name
        fields[n].idvar = ivar
        if printfield then print,'Field ',n,': ',fields[n].name
;
; Set dimension sizes and coordinate vars for lat,lon,lev,time:
;
        for idim = 0,var.ndims-1 do begin
          ncdf_diminq,cdfid,var.dim(idim),dimname,dimsize
          if (strpos(dimname,'time') gt -1) then begin
            fields[n].ntime=dimsize

;           print,'readfile: field ',fields[n].name,' ndims=',var.ndims,$
;             ' dimname=',dimname,' dimsize=',dimsize,' fields[n].ntime=',$
;             fields[n].ntime

            if (dimsize gt 0) then begin
              timeid = ncdf_varid(cdfid, dimname)
              times = fltarr(dimsize)
              ncdf_varget,cdfid,timeid,times
              if ptr_valid(fields[n].times) then $
                *fields[n].times=times else fields[n].times=ptr_new(times)
            endif
          endif
          if (strpos(dimname,'lat' ) gt -1) then begin
            fields(n).nlat=dimsize
            if (dimsize gt 0) then begin
              latid = ncdf_varid(cdfid, dimname)
              lats = fltarr(dimsize)
              ncdf_varget,cdfid,latid,lats
              if ptr_valid(fields[n].lats) then $
                *fields[n].lats=lats else fields[n].lats=ptr_new(lats)
            endif
          endif
          if (strpos(dimname,'lon' ) gt -1) then begin
            fields(n).nlon=dimsize
            if (dimsize gt 0) then begin
              lonid = ncdf_varid(cdfid, dimname)
              lons = fltarr(dimsize)
              ncdf_varget,cdfid,lonid,lons
              if ptr_valid(fields[n].lons) then $
                *fields[n].lons=lons else fields[n].lons=ptr_new(lons)
            endif
          endif
          if (strpos(dimname,'lev' ) gt -1) then begin
            fields(n).nlev=dimsize
            if (dimsize gt 0) then begin
              levid = ncdf_varid(cdfid, dimname)
              levs = fltarr(dimsize)
              ncdf_varget,cdfid,levid,levs
              if ptr_valid(fields[n].levs) then $
                *fields[n].levs=levs else fields[n].levs=ptr_new(levs)
;
; Get lev label (name and units) (assume lev var name is dimname dimension)
              long_name = ''
              units = ''
              for ivartmp = 0,file_inq.nvars-1 do begin
                vartmp = ncdf_varinq(cdfid,ivartmp)
                if vartmp.name eq dimname then begin
                  for ii=0,vartmp.natts-1 do begin
                    attname = ncdf_attname(cdfid,ivartmp,ii) 
                    if (attname eq 'long_name') then begin
                      ncdf_attget,cdfid,ivartmp,'long_name',rd_long_name
                      long_name = string(rd_long_name)
                    endif
                    if (attname eq 'units') then begin
                      ncdf_attget,cdfid,ivartmp,'units',rd_units
                      units = string(rd_units)
                    endif
                  endfor
                endif
              endfor
              fields[n].levname = long_name
              fields[n].levunits = units
              if strpos(long_name,"pressure") gt -1 or $
                 strpos(long_name,"PRESSURE") gt -1 or $
                 strpos(long_name,"Pressure") gt -1 then $
                 fields[n].levshortname = 'ZP'
              if strpos(long_name,"height") gt -1 or $
                 strpos(long_name,"Height") gt -1 then $
                 fields[n].levshortname = 'HEIGHT'
             ;print,'readfile: n=',n,' fields[n].levname=',fields[n].levname,$
             ;      'levunits=',fields[n].levunits
            endif
          endif
          fields[n].dimnames(idim) = dimname
        endfor
        if (var.natts gt 0) then begin
;
; Get field attributes:
;
          for i=0,var.natts-1 do begin
            attname = ncdf_attname(cdfid,ivar,i) 
            case attname of
            'long_name': begin
              ncdf_attget,cdfid,ivar,'long_name',long_name
              long_name = string(long_name)
              if strlen(strcompress(long_name,/remove_all)) eq 0 then $
                fields[n].long_name = fields[n].name $
              else fields[n].long_name = long_name
              if printfield then print,'  long_name: ',fields[n].long_name
            end
            'units': begin
              ncdf_attget,cdfid,ivar,'units',units
              fields[n].units = string(units)
              if printfield then print,'  units: ',fields[n].units
            end
            'differences': begin
              ncdf_attget,cdfid,ivar,'differences',difftype
              fields[n].difftype = string(difftype)
;             difftype = string(difftype)
;             if strpos(difftype,'PERC') ge 0 or strpos(difftype,'perc') ge 0 $
;               then fields[n].difftype = '%'
              if printfield then print,'  difftype: ',fields[n].difftype
            end
            'type': begin

            end
            else: begin
              print,'Did not read attribute "',attname,'" of field ',$
                fields[n].name
            end
            endcase
          endfor ; 0,var.natts-1
        endif ; var.natts > 0
;
; Set grid type (geographic or magnetic):
;
        fields[n].grid_type = 'geographic'
        for idim = 0,var.ndims-1 do begin
          if (strpos(fields[n].dimnames(idim),'mlon') gt -1) then $
            fields[n].grid_type = 'magnetic'
        endfor
;
; File from which field is read:
        fields[n].ncfile = info.file
      endif ; is lat,lon,time var [n]
    endif ; var.ndims > 0
  endfor ; ivar=0,file_inq.nvars-1
endif ; nvars > 0
return,fields
end
;-----------------------------------------------------------------------
pro readfile,info
;
; Read netcdf file attached to info.cdfid.
;
cdfid = info.cdfid
print,'Reading file ',info.file,'...'
;
; Check for difference fields:
;
diffs = 0
fileinfo = ncdf_inquire(cdfid)
for i=0,fileinfo.ngatts-1 do begin
  name = ncdf_attname(cdfid,/global,i)
  if name eq 'difference_fields' then begin
    ncdf_attget,cdfid,/global,name,diffbytes
;   print,'Difference fields = ',string(diffbytes)
    if string(diffbytes) eq 'yes' then begin
      diffs = 1
      print,'This file contains difference fields.'
    endif
  endif
endfor
info.diffs = diffs
;
; Get fields:
;
fields = getvars(info)
pfields = ptr_new(fields)
info.fields = pfields
info.nflds = n_elements(fields)
;
; Read model times:
;
idmtime = ncdf_varid(cdfid,'mtime')
if idmtime eq -1 then begin
  print,'>>> WARNING readfile: no mtime variable.'
  stop,'mtime'
endif
mtime_struct = ncdf_varinq(cdfid,idmtime)
id_ntimes = mtime_struct.dim[1] ; mtime[3,ntimes]
dum=''
ncdf_diminq,cdfid,id_ntimes,dum,ntimes
mtimes = intarr(3,ntimes)
ncdf_varget,cdfid,idmtime,mtimes
;
; Read missing data value (global attribute in netcdf file):
;
missing = ncdf_attinq(cdfid,/global,'missing_value')
if missing.length ne 1 then begin
  print,'>>> WARNING readfile: cannot find missing_value variable'
  print,'    Will default to missing_value = 1.e+36'
  missing_value = 1.e+36
endif else begin
  ncdf_attget,cdfid,/global,'missing_value',missing_value
; print,'Read missing_value = ',missing_value
endelse
;
; Check if this file was processed by f90 tgcmproc:
;
; A known IDL bug in 5.3 and 5.4 causes a crash on ncdf_attinq/attget:
; This appears to work with IDL 5.6.
;
; Opened file /vishnu/e/foster/tgcmproc/FOSTER.tgcm24r.sy2n005.nc
; % NCDF_ATTINQ: Attribute inquiry failed, name "processed" already in use.
; % NCDF_ATTGET: Attribute inquiry failed, name "processed" already in use.
; % NCDF_ATTGET: File format error. Attribute processed has a unknown type
; % Execution halted at:  READFILE          214 /home/foster/tgcmproc/readfile.pro
; %                       OPENFILE_EVENT    118 /home/foster/tgcmproc/tgcmproc.pro
; %                       WIDGET_PROCESS_EVENTS       
; %       
;
version = float(!version.release)
tgcmproc_file = 0 ; init
ncdf_attget,cdfid,/global,'model_version',model_version
model_version = string(model_version)
print,'model_version=',model_version
;
if version ge 5.6 then begin
  tgcmproc_att = ncdf_attinq(cdfid,/global,'tgcmproc')
  if tgcmproc_att.length eq 3 then begin
    ncdf_attget,cdfid,/global,'tgcmproc',value
    tgcmproc_file = 1
    print,'Note: This IS an f90 tgcmproc output file.'
  endif else begin
    tgcmproc_file = 0
    print,'Note: This is NOT an f90 tgcmproc output file.'
  endelse
endif else begin
;
; With this version (< 5.6) we cannot test for an attribute that
; does not exist (because of idl bug), so if Z is on the history,
; test its units: Z units will be CM on model files, and KM on 
; tgcmproc output files.
;
  zid = ncdf_varid(cdfid,'Z')
  if zid ge 0 then begin
    ncdf_attget,cdfid,zid,'units',zunits
    print,'zunits=',string(zunits)
    zunits = string(zunits)
    case zunits of
      'CM': begin
        tgcmproc_file = 0
        print,'Note: This is NOT an f90 tgcmproc output file ',$
          '(zunits=',string(zunits),')'
      end
      'KM': begin
        if model_version ne 'tgcm15' then begin
          tgcmproc_file = 1
          print,'Note: This IS an f90 tgcmproc output file ',$
          '(zunits=',string(zunits),')'
;
; tgcm15 files have erroneous KM units on Z:
        endif else begin
          tgcmproc_file = 0
          print,'Note: This is NOT an f90 tgcmproc output file ',$
          '(tgcm15 file with incorrect KM units for Z)'
        endelse
      end
      else: begin
        print,'>>> unknown zunits: ',zunits,' -- cannot determine',$
          ' if this file is a tgcmproc output file.'
        stop
      end
    endcase
  endif
;
; vnid = ncdf_varid(cdfid,'VN')
; if vnid ge 0 then begin
;   ncdf_attget,cdfid,vnid,'units',vn_units
;   print,'vn_units=',string(vn_units)
;   vn_units = string(vn_units)
;   case vn_units of
;     'CM/S': begin
;       tgcmproc_file = 0
;       print,'Note: This is NOT an f90 tgcmproc output file ',$
;         '(vn_units=',string(vn_units),')'
;     end
;     'M/S': begin
;       tgcmproc_file = 1
;       print,'Note: This IS an f90 tgcmproc output file ',$
;       '(vn_units=',string(vn_units),')'
;     end
;     else: begin
;       print,'>>> unknown vn_units: ',vn_units,' -- cannot determine',$
;         ' if this file is a tgcmproc output file.'
;       stop
;     end
;   endcase
; endif
endelse
;
; Get lons,lats,levs:
;
idlon = ncdf_varid(cdfid,'lon')
if idmtime gt -1 then ncdf_varget,cdfid,idlon,lons else lons = 0.
idlat = ncdf_varid(cdfid,'lat')
if idmtime gt -1 then ncdf_varget,cdfid,idlat,lats else lats = 0.
;print,'readfile: lons=',lons
;print,'readfile: lats=',lats
;
; Create file structure:
 file = {                  $
   path:info.file,         $ ; path to netcdf file
   cdfid:info.cdfid,       $ ; file id connected to file
   nflds:info.nflds,       $ ; number of 3d or 4d fields to display
   mtimes:mtimes,          $ ; model times (int(3))
   ntimes:ntimes,          $ ; number of model times
   diffs:diffs,            $ ; diffs==1 if difference fields (global file attribute)
   lons:lons,              $
   lats:lats,              $
   missing_value:missing_value, $ ; missing data value
   tgcmproc_file:tgcmproc_file  $
 }
pfile = ptr_new(file)
info.pfile = pfile
end