FUNCTION levToPlotLev, inLevs
;
;  Look at ln pressure levels above about -10 and determine the delta to use to form the standard
;  ln pressure grid to be used for plotting
;
  nLevsIn = N_ELEMENTS(inLevs)
  indsLow = WHERE(inLevs LE -10)
  IF indsLow[0] NE -1 THEN BEGIN

    nindsLow = N_ELEMENTS(indsLow)
    indLow = indsLow[nindsLow-1]
    levsLow = inLevs[indLow:nLevsIn-2]
    levsHi = inLevs[indLow+1:nLevsIn-1]
    delLevsPlot = MIN(levsHi - levsLow)
    nLevsPlot = FIX((inLevs[nLevsIn-1] - inLevs[0])/delLevsPlot) + 2

  ENDIF

  levLastPlot = inLevs[nLevsIn-1]
  levsPlot = FINDGEN(nLevsPlot)*delLevsPlot
  levsPlot = levsPlot - levsPlot[nLevsPlot-1]
  levsPlot = levsPlot + (levLastPlot - 0.001)
  
  RETURN, levsPlot

END
;
;--------------------------------------------------------------------------
;
FUNCTION levToLnPLev, inLevs
;
;  Flip the levels and convert to ln pressure
;
  nLevs = N_ELEMENTS(inlevs) 
  outLevs = inLevs

  FOR iLev=0,nLevs-1 DO BEGIN

    kLev = nLevs - iLev - 1
    outLevs(kLev) =  ALOG(5.0E-07/inlevs(iLev))  

  ENDFOR	      

  RETURN, outLevs

END
;
;--------------------------------------------------------------------------
;
FUNCTION getTimesW, cdFIDW
;
; Read WACCM times and get into TIMEGCM format:
;
  IDTimeD = ncdf_dimid(cdFIDW,'time')

  if IDTimeD eq -1 then begin

    print,'>>> WARNING readfile: no time dimension.'
    stop,'time'

  endif

  ncdf_diminq,cdFIDW,IDTimeD,timeName,nTimes

  times = INTARR(3,nTimes)
;
;  Need to get day, hour, minute for WACCM times
;
  IDDate = ncdf_varid(cdFIDW, 'date')
  ncdf_varget,cdFIDW,IDDate,inDates

  strDate = STRTRIM(STRING(inDates),2)
  inMonths = FIX(STRMID(strDate,4,2))
  inDays = FIX(STRMID(strDate,6,2))

  inYears = inDates/10000

  julDayBeg = JULDAY(1,1,inYears[0])

  julDayNow = JULDAY(inMonths,inDays,inYears)
  
  IDSec = ncdf_varid(cdFIDW, 'nscur')
  ncdf_varget,cdFIDW,IDSec,inSecs

  inHours = inSecs/3600.0

  inMinSecs = inHours - FIX(inHours)
 
  times[0,0:nTimes-1] = julDayNow - julDayBeg + 1
  times[1,0:nTimes-1] = FIX(inHours)
  times[2,0:nTimes-1] = FIX(inMinSecs*60)
;
;  Need to handle case of 1 time in file
;  
  if nTimes eq 1 then begin
  
    timestemp = intarr(3,2)
    timestemp(*,0) = times(*,0)
    timestemp(*,1) = -1
    times = timestemp
    
  endif

  RETURN, times

END
;
;--------------------------------------------------------------------------
function getvarsW,info
;
; Find variables on cdFIDW 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
cdFIDW = info.cdfid
ncdf_control,cdFIDW,/verbose
file_inq = ncdf_inquire(cdFIDW)
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(cdFIDW,ivar)
  if (var.ndims gt 0) then begin
    dimnames = ' '
    for idim = 0,var.ndims-1 do begin
      ncdf_diminq,cdFIDW,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,cdfidW,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(?)
  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)
  n = -1
  fields1D = replicate({field1dstruct},nvars1D)
  n1D = -1
  levNotFound = 0
  for ivar = 0,file_inq.nvars-1 do begin
    var = ncdf_varinq(cdfidW,ivar)
    if (var.ndims gt 0) then begin
      dimnames = ' '
      for idim = 0,var.ndims-1 do begin
        ncdf_diminq,cdfidW,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 is not the first field and no field with vertical levels has been found and this field
;  has vertical levels, 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
;
; Set dimension sizes and coordinate vars for lat,lon,lev,time:
;
        for idim = 0,var.ndims-1 do begin
          ncdf_diminq,cdfidW,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
	      
  	      if info.nfiles gt 1 then begin
	      
                ntimes = 0
	        files = *info.files
	        cdfid = cdFIDW
	      
		for ifile = 0, info.nfiles-1 do begin
     
                  if ifile gt 0 then cdFID = ncdf_open(files[ifile],/nowrite)
                  timeid = ncdf_dimid(cdfid,'time')
                  ncdf_diminq,cdFID,timeid,name,ntime
                  ncdf_varget,cdFID,name,timesIn
		  timesAllIn[ntimes:ntimes+ntime-1] = timesIn[*]
                  ntimes = ntimes + ntime
                  if ifile gt 0 then ncdf_close,cdFID
 		
		endfor

                fields[n].ntime=ntimes
                times = timesAllIn

                cdfidW = info.cdfid ;Set file id back to first file
		
	      endif else begin	         

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

              endelse
 
              model_name = 'WACCM'

              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(cdfidW, dimname)
              lats = fltarr(dimsize)
              ncdf_varget,cdfidW,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(cdfidW, dimname)
              lons = fltarr(dimsize)
              ncdf_varget,cdfidW,lonid,lons

              convertlons,info,lons
;
;  Lons are now in TIME-GCM order but need the values adjusted
;
              iHighLons = WHERE(lons ge 180.0)
              lons(iHighLons) = lons(iHighLons) - 360.0
              
              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(cdfidW, dimname)
              levs = fltarr(dimsize)
              ncdf_varget,cdfidW,levid,rdLevs
;
;  Flip vertical axis to go from ground up and convert to ln pressure
;
              levs = levToLnPLev(rdLevs)

              if ptr_valid(fields[n].levs) then $
                *fields[n].levs=levs else fields[n].levs=ptr_new(levs)
;
;  Form a regular ln pressure grid to use for plotting and save to field structure
;
              levsPlot = levToPlotLev(levs)

              if ptr_valid(fields[n].levsPlot) then $
                *fields[n].levsPlot = levsPlot else fields[n].levsPlot = ptr_new(levsPlot)
;
; 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(cdfidW,ivartmp)
                if vartmp.name eq dimname then begin
                  for ii=0,vartmp.natts-1 do begin
                    attname = ncdf_attname(cdfidW,ivartmp,ii) 
                    if (attname eq 'long_name') then begin
                      ncdf_attget,cdfidW,ivartmp,'long_name',rd_long_name
                      long_name = string(rd_long_name)
                    endif
                    if (attname eq 'units') then begin
                      ncdf_attget,cdfidW,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 ;lev if
          fields[n].dimnames(idim) = dimname
	  
	  if (strpos(dimname,'ilev' ) gt -1) then fields[n].levstype = 'interface'
	  
        endfor ;dims loop
        if (var.natts gt 0) then begin
;
; Get field attributes:
;
          for i=0,var.natts-1 do begin
            attname = ncdf_attname(cdfidW,ivar,i) 
            case attname of
            'long_name': begin
              ncdf_attget,cdfidW,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,cdfidW,ivar,'units',units
              fields[n].units = string(units)
              if printfield then print,'  units: ',fields[n].units
            end
            'differences': begin
              ncdf_attget,cdfidW,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,cdfidW,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,cdfidW,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 = cdfidW
	
	  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]

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

        endif else begin ; Only one input file	   

          fields1D[n1D].ntime=dimsize

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

        endelse

        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(cdfidW,ivar,i) 
            case attname of
            'long_name': begin
              ncdf_attget,cdfidW,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,cdfidW,ivar,'units',units
              fields1D[n1D].units = string(units)
              if printfield then print,'  units: ',fields1D[n1D].units
            end
            'differences': begin
              ncdf_attget,cdfidW,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,cdfidW,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
;-----------------------------------------------------------------------
pro readFileW,info
;
; Read WACCM netcdf history file attached to info.cdfid.
;
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 

cdFIDW = info.cdfid
print,'Reading file ',info.file,'...'
;
; Check for difference fields:
;
diffs = 0
fileinfo = ncdf_inquire(cdFIDW)

;help,fileinfo,/struct
;exit

for i=0,fileinfo.ngatts-1 do begin
  name = ncdf_attname(cdFIDW,/global,i)
  if name eq 'difference_fields' then begin
    ncdf_attget,cdFIDW,/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 done for TIMEGCM so keep variable here:
;
oldhist = 0
tgcmproc_file = 0
;
; Get fields:
;
fields = getvarsW(info)
pfields = ptr_new(fields)
info.fields = pfields
info.nflds = n_elements(fields)
;
;  Open file(s) and get number 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 cdFIDW = ncdf_open(files[ifile],/nowrite)

    id = ncdf_dimid(cdFIDW,'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,cdFIDW,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

    mtime = getTimesW(cdFIDW)
;
;  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,cdFIDW

  endfor ; ifile=0,nfiles-1

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

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

endif else begin ; One input file 
 
  times = getTimesW(cdFIDW)
  timeInfo = SIZE(times)
  nTimes = timeInfo(2)
  if nTimes eq 2 and times[0,nTimes-1] eq -1 then nTimes = 1
  mtimes = times[*,0:nTimes-1]
  mtime_files[*,0:ntimes-1,0] = times[*,0:nTimes-1]
  ntime_files = nTimes
  mtime_beg[*] = mtimes[*,0]
  ifile_beg = 0
  mtime_end[*] = mtimes[*,ntimes-1]
  ifile_end = 0

endelse
;
; Get lons,lats,levs:
;
IDLon = ncdf_varid(cdFIDW,'lon')
if IDLon gt -1 then BEGIN

  ncdf_varget,cdFIDW,IDLon,lons

  convertlons, info, lons
;
;  Lons are now in TIME-GCM order but need the values adjusted
;
  iHighLons = WHERE(lons ge 180.0)
  lons(iHighLons) = lons(iHighLons) - 360.0

ENDIF else BEGIN 

  lons = 0.

ENDELSE

IDLat = ncdf_varid(cdFIDW,'lat')
if IDLat gt -1 then ncdf_varget,cdFIDW,IDLat,lats else lats = 0.
;print,'readfile: lons=',lons
;print,'readfile: lats=',lats

IDLev = ncdf_varid(cdFIDW,'lev')
if IDLev gt -1 then BEGIN 

  ncdf_varget,cdFIDW,IDLev,rdLevs
;print,'readfile: rdlevs =',rdlevs

  levs = levToLnPLev(rdLevs)

ENDIF else BEGIN 

  levs = 0.

ENDELSE

IDLev = ncdf_varid(cdFIDW,'ilev')
if IDLev gt -1 then BEGIN 

  ncdf_varget,cdFIDW,IDLev,rdILevs

  iLevs = levToLnPLev(rdILevs)

ENDIF else BEGIN

  ilevs = 0.

ENDELSE
;print,'readfile: levs (midpoints)=',levs
;print,'readfile: ilevs (interfaces)=',ilevs
;
missing_value = 1.0E+36

; 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