;
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
nvars1D = 0  ;Counter for 1D fields in file
maxtimes = 10000            ; max number of total times (all files)
if info.nfiles gt 1 then timesAllIn  = FLTARR(maxtimes)  ; times from all input files

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
;  
;  Count 1D time fields in file
;   
  IF var.ndims EQ 1 THEN BEGIN
    ncdf_diminq,cdfid,var.dim[0],dimname,dimsize
    IF dimname EQ 'time' THEN BEGIN
      nvars1D = nvars1D + 1
    ENDIF     
  ENDIF 
   
endfor ; ivar=0,nvars-1
;
; Define a structure to hold 1D (time only for now) variables from the file:
;
; Structure for a variable with time dimension.
;
field1D = {field1dstruct,            $ ; a 1D field structure type
  name:' ',                          $ ; short name (char8)
  long_name: ' ',                    $ ; long name (may be blank)
  units: ' ',                        $ ; units (may be blank)
  tunits: ' ',                       $ ; time units
  idvar:0,                           $ ; netcdf variable id
  ntime:0,                           $ ; time dimension size
  dimnames:strarr(1),                $ ; dimension name in
  ncfile:'',                         $ ; file from which field is read
  difftype:'',                       $ ; difference type ('' or '%')
  missing_value:0.,                  $ ; missing data value for this field
  times:ptr_new(/allocate_heap),     $ ; times coord var
  data:ptr_new(/allocate_heap)       $ ; the data
}
;
; 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)
  nlevlp:0,                          $ ; vertical dimension size for ln pressure (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)
  levstype:'midpoint',               $ ; levels coord type, either on midpoints or interfaces
  levshortname:'',                   $ ; short name of lev coord var
  levsPlot:ptr_new(/allocate_heap),  $ ; ln pressure level regular grid for plotting(Consistent with WACCM)
  missing_value:0.,                  $ ; missing data value for this field
  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)
  print, 'Making fields structure size ', nvars
  n = -1      
  fields1D = replicate({field1dstruct},nvars1D)
  n1D = -1
  levNotFound = 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
      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
;
;  Need first field to be 4D including vertical levels so check and if not the case switch
;  with a later field which has vertical levels
;
;  See if this field has vertical levels in dimension list
;
        lev_substr = strpos(dimnames,'lev')
;
;  If first field and doesn't have vertical levels then put in second element of field structure
;
	IF ( n EQ 0 AND lev_substr LT 0 ) THEN BEGIN
	
	  n = n + 1
	  levNotFound = 1
;
;  If this field has vertical levels, this is not the first field, and no field with vertical 
;  levels has been found then put this field in first element of field structure
;	  
	ENDIF ELSE IF ( n GT 0 AND levNotFound EQ 1 AND lev_substr GE 0 ) THEN BEGIN
	
	  nSave = n
	  n = 0
	  levNotFound = 2
;
;  If first field with levels was found and put in first element of field array then continue filling
;  field structure
;	  
	ENDIF ELSE IF ( n EQ 1 AND levNotFound EQ 2 ) THEN BEGIN
	
	  n = nSave
	  
	ENDIF
;
; Set name and netcdf var id:
;
        fields[n].name = var.name
        fields[n].idvar = ivar
        if printfield then print,'Field ',n,': ',fields[n].name
        ncdf_attget,cdfid,/global,'model_name',model_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
	    
;           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

  	      if info.nfiles gt 1 then begin
	      
	        ntimes = 0
                files = *info.files
	        ncdfid = cdfid
	      
		for ifile = 0, info.nfiles-1 do begin
     
                  if ifile gt 0 then ncdfid = ncdf_open(files[ifile],/nowrite)
                  timeid = ncdf_dimid(ncdfid,'time')
                  ncdf_diminq,ncdfid,timeid,name,ntime
                  ncdf_varget,ncdfid,name,timesIn
		  timesAllIn[ntimes:ntimes+ntime-1] = timesIn[*]
                  ntimes = ntimes + ntime
 		  if ifile gt 0 then ncdf_close,ncdfid
		 
		endfor
;
; Save number of times in field structure and point back to first input file id
;
                fields[n].ntime=ntimes
                times = timesAllIn
		
	      endif else begin	         
                
                fields[n].ntime=dimsize
                timeid = ncdf_varid(cdfid, dimname)
                times = fltarr(dimsize)
		ncdf_varget,cdfid,timeid,times

              endelse

;creates the rotation of time data when using a mars model
              
;             ncdf_attget,cdfid,/global,'model_name',model_name
              model_name=string(model_name)
              if model_name eq 'mtgcm' then begin
                    times = times + (60*12)
              endif 
              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 or $
                 strpos(long_name,"Pressure") gt -1 then $
                 fields[n].levshortname = 'ZP'
;
; New history file format (12/05) distinguishes between midpoint
; and interface levels:
              if strpos(long_name,"midpoint") gt -1 or $
                 strpos(long_name,"interface") gt -1 then $
                 fields[n].levshortname = 'ZP'
;
; Levels will be height only if file is tgcmproc_f90 post-proc output file:
              if strpos(long_name,"height") gt -1 or $
                 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
	  if (strpos(dimname,'ilev' ) gt -1) then fields[n].levstype = 'interface'
        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
;
; Only new histories have missing_value on per field basis:
            'missing_value': begin
              ncdf_attget,cdfid,ivar,'missing_value',missing_value
              fields[n].missing_value = missing_value
              if printfield then print,'  missing_value: ',fields[n].missing_value
            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
    
    IF var.ndims EQ 1 THEN BEGIN
   
      ncdf_diminq,cdfid,var.dim[0],dimname,dimsize

      IF dimname EQ 'time' THEN BEGIN

        n1D = n1D + 1      

; Set name and netcdf var id:
;
        fields1D[n1D].name = var.name
        fields1D[n1D].idvar = ivar
        if printfield then print,'Field ',n1D,': ',fields1D[n1D].name

  	if info.nfiles gt 1 then begin
	
	  ntimes = 0
          files = *info.files
	  ncdfid = cdfid
	
	  for ifile = 0, info.nfiles-1 do begin
     
            if ifile gt 0 then ncdfid = ncdf_open(files[ifile],/nowrite)
            timeid = ncdf_dimid(ncdfid,'time')
            ncdf_diminq,ncdfid,timeid,name,ntime
            ncdf_varget,ncdfid,name,timesIn
	    timesAllIn[ntimes:ntimes+ntime-1] = timesIn[*]
            ntimes = ntimes + ntime
 	    if ifile gt 0 then ncdf_close,ncdfid
	   
	  endfor
;
; Save number of times in field structure and point back to first input file id
;
          fields1D[n1D].ntime=ntimes
          times = timesAllIn[0:ntimes-1]

          cdfid = info.cdfid ;Set file id back to first file

        endif else begin	   

          fields1D[n1D].ntime=dimsize

          timeid = ncdf_varid(cdfid, dimname)
          times = fltarr(dimsize)
          ncdf_varget,cdfid,timeid,times

        endelse

;creates the rotation of time data when using a mars model
              
        ncdf_attget,cdfid,/global,'model_name',model_name
        model_name=string(model_name)
        if model_name eq 'mtgcm' then begin
          times = times + (60*12)
        endif 
        if ptr_valid(fields1D[n1D].times) then $
                    *fields1D[n1D].times=times else fields1D[n1D].times=ptr_new(times)

        fields1D[n1D].dimnames[0] = dimname

        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 $
                fields1D[n1D].long_name = fields1D[n].name $
              else fields1D[n1D].long_name = long_name
              if printfield then print,'  long_name: ',fields1D[n1D].long_name
            end
            'units': begin
              ncdf_attget,cdfid,ivar,'units',units
              fields1D[n1D].units = string(units)
              if printfield then print,'  units: ',fields1D[n1D].units
            end
            'differences': begin
              ncdf_attget,cdfid,ivar,'differences',difftype
              fields1D[n1D].difftype = string(difftype)
;             difftype = string(difftype)
;             if strpos(difftype,'PERC') ge 0 or strpos(difftype,'perc') ge 0 $
;               then fields1D[n1D].difftype = '%'
              if printfield then print,'  difftype: ',fields1D[n1D].difftype
            end
            'type': begin

            end
;
; Only new histories have missing_value on per field basis:
            'missing_value': begin
              ncdf_attget,cdfid,ivar,'missing_value',missing_value
              fields1D[n1D].missing_value = missing_value
              if printfield then print,'  missing_value: ',fields1D[n1D].missing_value
            end
            else: begin
;             print,'Did not read attribute "',attname,'" of field ',$
;               fields1D[n1D].name
            end
            endcase
          endfor ; 0,var.natts-1
        endif ; var.natts > 0

        IF fields1D[n1D].name EQ 'time' THEN fields1D.tunits = fields1D[n1D].units

;
; File from which field is read:
        fields1D[n1D].ncfile = info.file
     
      ENDIF  ;Only time dimension     
    ENDIF  ;1D variable
    
  endfor ; ivar=0,file_inq.nvars-1
  
  pfields1D = ptr_new(fields1D)
  info.fields1D = pfields1D  
  info.nflds1D = n_elements(fields1D)
  
; if model_name eq 'mtgcm' then $
;   print,'read mars global attributes adding rotation to minutes'
endif ; nvars > 0
return,fields
end ; function getvars
;-----------------------------------------------------------------------
pro readfile,info

maxtimes = 10000            ; max number of total times (all files)
mtimes  = FLTARR(3,maxtimes)  ; times from all input files
maxtime_file = 500	  ; max number of times per file
mtime_files = intarr(3,maxtime_file,info.nfiles)  ; mtimes on each file
if info.nfiles eq 1 then mtime_files = intarr(3,maxtime_file,info.nfiles+1)  ; mtimes on each file
ntime_files = intarr(info.nfiles)		    ; number of times on each file
ntimes = 0
mtime_beg = intarr(3)  ; beginning input model time 
mtime_end = intarr(3)  ; ending input model time 
;
; 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)

;help,fileinfo,/struct
;exit

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
;
; 10/28/05: Check for old history:
;
oldhist = 1
for i=0,fileinfo.nvars-1 do begin
  var = ncdf_varinq(cdfid,i)
  if (var.name eq 'LBC') then oldhist = 0 ; is new hist
endfor
;print,'oldhist = ',oldhist
;
; Get fields:
;
fields = getvars(info)
 
pfields = ptr_new(fields)
info.fields = pfields
info.nflds = n_elements(fields)
;
; Read model times:
;
; Open file(s) and get number of times (histories):
;
if info.nfiles gt 1 then begin

  files = *info.files
  mtimestemp = mtimes

  for ifile=0,info.nfiles-1 do begin

    if ifile gt 0 then cdfid = ncdf_open(files[ifile],/nowrite)

    id = ncdf_dimid(cdfid,'time')
    if (id eq -1) then begin
      print,'>>> Open_files: Could not find id of time dimension on file ',files[ifile]
      return
    endif
    name = ''
    ncdf_diminq,cdfid,id,name,ntime
    if (ntimes+ntime gt maxtimes) then begin
      print,'>>> Open_files: ntimes gt maxtimes: ntimes=',ntimes,' maxtimes=',maxtimes
      return
    endif 
    if (ntime gt maxtime_file) then begin
      print,'>>> Open_files: ntime gt maxtime_file: ntime=',ntime,' maxtime_file=',$
        maxtime_file
      return
    endif
;
; Get model times and add to total mtimes:
;
    id = ncdf_varid(cdfid,'mtime')
    if (id eq -1) then begin
      print,'>>> Open_files: Could not find id of mtime var on file ',files[ifile]
      return
    endif
  
    ncdf_varget,cdfid,id,mtime
;
;  Get first and last time and first and last file index of first and last input files
;    
    if ifile eq 0 then begin
      mtime_beg[*] = mtime[*,0]
      ifile_beg = ifile
    endif
    if ifile eq info.nfiles-1 then begin
      mtime_end[*] = mtime[*,ntime-1]
      ifile_end = ifile
    endif
    
    mtimestemp[*,ntimes:ntimes+ntime-1] = mtime[*,*]
    ntimes = ntimes + ntime

    if (ntime gt 1) then begin
      for it=1,ntime-1 do begin
        deltamin = mtime_to_mins(mtime[*,it])-mtime_to_mins(mtime[*,it-1])
        if deltamin ne 0 then break
      endfor
    endif else deltamin = 0.
    mtime_files[*,0:ntime-1,ifile] = mtime[*,*]
    ntime_files[ifile] = ntime
    print,format="('File ',a,' ntime=',i4,' (',3(i3,','),' to ',3(i3,','),' by ',i4,')')",$
      file_basename(files[ifile]),ntime,mtime[*,0],mtime[*,ntime-1],deltamin
    if ifile gt 0 then ncdf_close,cdfid

  endfor ; ifile=0,nfiles-1

  mtimes = mtimestemp[*,0:ntimes-1]

  cdfid = info.cdfid ;Set file id back to first file

endif else begin ; One input file 
  
  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
  mtime_files[*,0:ntimes-1,0] = mtimes[*,0:ntimes-1]
  ntime_files = nTimes
  mtime_beg[*] = mtimes[*,0]
  ifile_beg = 0
  mtime_end[*] = mtimes[*,ntimes-1]
  ifile_end = 0

endelse
;
; 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
endelse

;prefill = ncdf_attinq(cdfid,/global,'_FillValue')
;print,'Read _FillValue = prefill = ',prefill

;
; 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 idlon gt -1 then ncdf_varget,cdfid,idlon,lons else lons = 0.
idlat = ncdf_varid(cdfid,'lat')
if idlat gt -1 then ncdf_varget,cdfid,idlat,lats else lats = 0.
;print,'readfile: lons=',lons
;print,'readfile: lats=',lats

idlev = ncdf_varid(cdfid,'lev')
if idlev gt -1 then ncdf_varget,cdfid,idlev,levs else levs = 0.
idlev = ncdf_varid(cdfid,'ilev')
if (idlev gt -1) then ncdf_varget,cdfid,idlev,ilevs else ilevs = 0.
;print,'readfile: levs (midpoints)=',levs
;print,'readfile: ilevs (interfaces)=',ilevs

;corrects data for mars .nc by shifting mtime by 12 hrs
ncdf_attget,cdfid,/global,'model_name',model_name
model_name=string(model_name)
if model_name eq 'mtgcm' then begin
; mtimes[1,*]=mtimes[1,*]+12
; dims = size(mtimes)
; if dims[0] eq 1 then dims=dims[0]
; if dims[0] eq 2 then dims = dims[2]
; for i=0,dims-1 do begin 
;   if mtimes[1,i] ge 24 then begin
;     mtimes[0,i] = mtimes[0,i] + 1
;     mtimes[1,i]= mtimes[1,i] -24
;   endif
; endfor
; print,'read mars global attributes adding rotation to model time'
endif
;
; 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,              $
   levs:levs,              $ ; midpoint levels
   ilevs:ilevs,            $ ; interface levels
   missing_value:missing_value, $ ; missing data value
   tgcmproc_file:tgcmproc_file, $
   oldhist:oldhist         $
}
 
pfile = ptr_new(file)
info.pfile = pfile
;
; Update state:
;
info.ntimes = ntimes
*info.mtimes = mtimes[*,0:ntimes-1]
*info.mtime_files = mtime_files
*info.ntime_files = ntime_files
info.mtime_beg = mtime_beg
info.ifile_beg = ifile_beg
info.mtime_end = mtime_end
info.ifile_end = ifile_end

end