;
; User input (called from main pro utproc below):
;
pro input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,$
  zpbot,zptop,mtime_beg,mtime_end
;
; Directories and filename roots of nc history files to read.
; (dirs and fileroots must have same number of elements)
;
 dirs = ["/oreo/d/foster/ganglu"]
 fileroots = ["GANGLU.tiegcm1.8.s_nov04_"] ; use form="(i03)"
;
; Beginning and ending file numbers to read:
;
 file_beg = [21]  & file_end = [25] ; ganglu tiegcm1.8.s_nov04
;
; Format of file numbers (2 or 3 digits):
;
 form = "(i02)"
;form = "(i03)"
;
; Beginning and ending model times:
;   If mtime_beg[0] day < 0 then use the first time on first file
;   If mtime_end[0] day < 0 then use the last time on last file
;
 mtime_beg = [314,0,10] & mtime_end  = [-1,0,0]
;
; 4d fields to read and plot:
;
 fields = ["TN","UN","VN","PEDERSEN","HALL","QWIND","QAMIE","FWIND","FAMIE","TEC","QJOULE","QJH_TOT"]
;
; Lons and lats of locations to read (loc_lons and loc_lats must
; have same number of elements). Nearest gridpoints will be chosen.
;
;loc_lons = [-180.,-180.,-180.,-180.,-180.,-180.,-180.,-180.,-180.]
 loc_lons = [0.,0.,0.,0.,0.,0.,0.,0.,0.]
 loc_lats = [  -80. , -70., -60., -50., 0., 50.,  60.,  70.,  80.]
;
; zpbot,zptop: Bottom, top zp pressures for y-axis (bottom < top):
; Can set zpbot=zptop to default to full vertical grid.
; If zpbot < bottom grid point, then bottom grid point will be used
; If zptop > top grid point, then top grid point will be used
;
 zpbot = -17. & zptop = 5.
;zpbot = -17. & zptop = -7.
;zpbot = -999. & zptop = 999.
;zpbot = -999. & zptop = 5.
end ; pro input
;
;-----------------------------------------------------------------------
; 9/07: This clrbar is adapted from /home/tgcm/tgcmproc_idl/clrbar.pro
;
pro clrbar,bw,bh,xo,yo,labels,BOTCOLOR=botcolor,TOPCOLOR=topcolor,$
  LABCOLOR=labcolor,CHARSIZE=charsize,LABLEFT=lableft,LABTOP=labtop,$
  LABOFFSET=laboffset,COLORTAB=colortab
;
if (bw le 0. or bw gt 1. or bh le 0. or bh gt 1. or $
    xo lt 0. or xo gt 1. or yo lt 0. or yo gt 1.) then begin
  print,'>>> clrbar: bad arg (all coords and sizes ',$
    'must be in normalized coords (0->1.)): bw,bh=',bw,bh,$
    ' xo,yo=',xo,yo
  return
endif
if (keyword_set(charsize)) then chsize=charsize else chsize = 1.5
if (keyword_set(labcolor)) then iclab = labcolor else iclab = !p.color
if (keyword_set(topcolor)) then itopclr = topcolor else $
  itopclr = !d.table_size-1
if (keyword_set(botcolor)) then ibotclr = botcolor else ibotclr = 0
nc = itopclr-ibotclr+1
if (nc le 0) then begin
  print,'>>> clrbar: no colors?? nc=',nc,' bot,top color indices=',$
    ibotclr,itopclr
  return
endif
nlab = n_elements(labels)
if (nlab gt 0 and nlab ne nc+1) then begin
  print,' '
  print,'>>> clrbar: please provide a label for each color plus 1'
  print,'    (labels can be empty strings): ncolors = ',nc
  return
endif
if (keyword_set(laboffset)) then begin
  laboff = laboffset
endif else begin
  laboff = .035
  if (keyword_set(labtop)) then laboff = .025
endelse
;
; Horizontal color bar:
;
if (bw ge bh) then begin
  cellw = bw/float(nc)
  ii = 0
  for i=ibotclr,itopclr do begin
    x0 = xo + float(ii)*cellw
    x1 = x0+cellw
    x = [x0,x0,x1,x1]
    y = [yo,yo+bh,yo+bh,yo]
    polyfill,x,y,color=i,/norm
    ii = ii+1
  endfor
;
; Labels:
;
  ii = 0
  for i=ibotclr,itopclr+1 do begin
    x0 = xo + float(ii)*cellw
    if (nlab gt 0) then begin
      if (keyword_set(labtop)) then yl = yo+bh+laboff else yl = yo-laboff
      xyouts,x0,yl,labels(ii),/norm,align=.5,charsize=chsize
      plots,[x0,x0],[yo,yo+bh],/norm,color=iclab
    endif
    ii = ii+1
  endfor
endif else $
  print,'>>> this code makes horizontal color bar only:',$
   ' please make bw > bh'
;
; Border:
;
  x = [xo,xo,xo+bw,xo+bw,xo]
  y = [yo,yo+bh,yo+bh,yo,yo]
  plots,x,y,color=iclab,/norm
end
;-----------------------------------------------------------------------
function ixfind, z,val,del
;
; Find nearest index to val in z, with increment del
; If not found, return -1
;
nz = n_elements(z)
if del gt 0 then begin
  for i=0,nz-1 do begin
    if (val ge z(i)-0.5*del and val le z(i)+0.5*del) then return,i
  endfor
endif else begin
  for i=0,nz-1 do begin
    if (val eq z(i)) then return,i
  endfor
endelse
return,-1
end
;-----------------------------------------------------------------------
pro dayofyear_to_monthday, dayofyear,imonth,iday
;
; Given day of year (1-365), return month and day
;
iday = 0
imonth = 0
if (dayofyear eq 0) then dayofyear = 1 ; to allow for model day of 0
if (dayofyear lt 1 or dayofyear gt 365) then begin
  print,'>>> bad dayofyear=',dayofyear
  return
endif
;                J, F, M, A, M, J, J, A, S, O, N, D
dayspermonth = [31,28,31,30,31,30,31,31,30,31,30,31]
nmonths = n_elements(dayspermonth)
;
id = 0
for m=0,nmonths-1 do begin
  for i=1,dayspermonth[m] do begin
    id = id+1
    if (id eq dayofyear) then begin
      imonth = m+1
      iday = i
      break  
    endif
  endfor
  if (imonth ne 0) then break ; day was found above
endfor
;print,'dayofyear_to_monthday: dayofyear=',dayofyear,' imonth=',imonth,' iday=',iday
end
;-----------------------------------------------------------------------
function monthday_to_dayofyear, imonth,iday
;                J, F, M, A, M, J, J, A, S, O, N, D
dayspermonth = [31,28,31,30,31,30,31,31,30,31,30,31]
day = 0
for i=0,11 do begin
  if (i+1 eq imonth) then begin
    day = day+iday
    break
  endif
  day = day+dayspermonth[i]
endfor
;print,'monthday_to_dayofyear: imonth=',imonth,' iday=',iday,' dayofyear=',day
return,day
end
;-----------------------------------------------------------------------
function mtime_to_mins,mtime ; model time to integer minutes
return,mtime(0)*24*60+mtime(1)*60+mtime(2)
end
;-----------------------------------------------------------------------
function mtime_to_ddays,mtime ; model time to decimal days
return,float(mtime(0))+float(mtime(1))/24.+float(mtime(2))/(24.*60.)
end
;-----------------------------------------------------------------------
function calcslt,ut,glon
  slt = ut+glon/15.
  if (slt lt 0.) then slt = slt+24.
  if (slt ge 24.) then slt = slt-24.
  return,slt
end
;-----------------------------------------------------------------------
function timeaxis,years,mtimes,deltamin,xtitle,nxticks
;
; Get julian day of first time:
;
dayofyear_to_monthday,mtimes(0,0),imonth0,iday0
julian0 = julday(imonth0,iday0,years[0],mtimes(1,0),mtimes(2,0),0)
;
; Generate x-axis coords in julian days (step size in minutes):
; (this array is returned as the function value)
;
ntimes = n_elements(years)
time = timegen(ntimes,start=julian0,step_size=deltamin,units="Minutes")
;
; Add 0.5 to times since the julian day starts at 12 pm noon 
; (s.a., caldat call below)
;
time = time+0.5 
;
; Construct x-axis label giving beginning,ending model times,
; delta time in miniutes, and total time in days:
;
ddays0 = mtime_to_ddays(mtimes[*,0])
ddays1 = mtime_to_ddays(mtimes[*,ntimes-1])
total_ddays = ddays1-ddays0 ; total number of decimal days 
nxticks = 0 ; let idl select
;if (total_ddays ge 364) then nxticks = 7 ; does not appear to work
xtitle = string(format=$
  "('Model time (yyyy:ddd:hh:mm) ',i4,':',i3.3,':',i2.2,':',i2.2,' to ')",$
  years[0],mtimes[*,0])
xtitle = xtitle + string(format=$
  "(i4,':',i3.3,':',i2.2,':',i2.2,' by ',i6,' minutes ')",$
  years[ntimes-1],mtimes[*,ntimes-1],deltamin)
xtitle = xtitle + string(format="('(',f7.2,' days)')",total_ddays)
return,time
end
;-----------------------------------------------------------------------
; Main:
;
pro utproc,ps=ps
input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,$
  zpbot,zptop,mtime_beg,mtime_end
ndirs = n_elements(dirs)
nfiles = 0
for idir=0,ndirs-1 do begin
  nfiles_dir = file_end[idir]-file_beg[idir]+1
  nfiles = nfiles+nfiles_dir
endfor
files = strarr(nfiles)
ifiles = 0
for idir=0,ndirs-1 do begin
  nfiles_dir = file_end[idir]-file_beg[idir]+1
  for ifile=0,nfiles_dir-1 do begin
    i = file_beg[idir]+ifile
    strnum = string(format=form,i)
    if (idir eq 0 and ifile eq 0) then file0 = fileroots[idir]+strnum+".nc"
    if (idir eq ndirs-1 and ifile eq nfiles_dir-1) then $
      file1 = fileroots[idir]+strnum+".nc"
    files[ifiles] = dirs[idir] + "/" + fileroots[idir] + strnum + ".nc"
    ifiles = ifiles+1 
  endfor
endfor
;
nfields = n_elements(fields)
f_units = strarr(nfields)
ncids = intarr(nfiles)
if (n_elements(loc_lons) ne n_elements(loc_lats)) then begin
  print,'>>> Number of loc_lons must equal number of loc_lats:'
  print,'    nloc_lons=',n_elements(loc_lons),' nloc_lats=',$
    n_elements(loc_lats)
  return
endif
nlocs = n_elements(loc_lons)
print,format="('nfiles=',i4,' nfields=',i4,' nlocs=',i4)",$
  nfiles,nfields,nlocs
;
mtime0 = mtime_beg
mtime1 = mtime_end
;
; Search for start and end times:
;
ifile0 = -1 & ifile1 = -1
ihist0 = -1 & ihist1 = -1
for ifile=0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  id = ncdf_varid(ncids[ifile],'mtime')
  if (id eq -1) then begin
    print,'>>> could not find id of mtime var on file ',files[ifile]
    return
  endif
  ncdf_varget,ncids[ifile],id,mtime
  ntime = n_elements(mtime)/3
  for i=0,ntime-1 do begin
;   print,format=$
;     "('Searching for mtime0=',3i4,' or mtime1=',3i4,' Found ',3i4,' on file ',a)",$
;     mtime0,mtime1,mtime[*,i],file_basename(files[ifile])
;
; If start day is < 0, use first history on first file:
    if (mtime0[0] lt 0 and ifile eq 0 and i eq 0) then begin
      ifile0 = ifile
      ihist0 = i
      print,format="('Using start mtime ',3i4,' (first history of first file ',a,')')",$
        mtime[*,ihist0],file_basename(files[ifile0])
    endif else begin
      if (mtime[0,i] eq mtime0[0] and mtime[1,i] eq mtime0[1] and $
          mtime[2,i] eq mtime0[2]) then begin
        ifile0 = ifile
        ihist0 = i
        print,format="('Found start mtime ',3i4,' on file ',a,' (history ',i3,')')",$
          mtime0,file_basename(files[ifile]),ihist0
      endif
    endelse
;
; If end day is < 0 or > 366, use last history of last file:
    if ((mtime1[0] lt 0 or mtime1[0] gt 366) and $
         ifile eq nfiles-1 and i eq ntime-1) then begin
      ifile1 = nfiles-1
      ihist1 = ntime-1
      print,format="('Using ending mtime ',3i4,' (last history on last file ',a,')')",$
        mtime[*,ihist1],file_basename(files[ifile1])
    endif else begin
      if (mtime[0,i] eq mtime1[0] and mtime[1,i] eq mtime1[1] and $
          mtime[2,i] eq mtime1[2]) then begin
        ifile1 = ifile
        ihist1 = i
        print,format="('Found end   mtime ',3i4,' on file ',a,' (history ',i3,')')",$
          mtime1,file_basename(files[ifile]),ihist1
      endif
    endelse
    if (ifile0 gt -1 and ifile1 gt -1) then break
  endfor ; i=0,ntime-1
  ncdf_close,ncids[ifile]
  if (ifile0 gt -1 and ifile1 gt -1) then break
endfor ; mtimes search: ifile=0,nfiles-1
if (ifile0 eq -1) then begin
  print,format="('>>> Could not find start mtime ',3i4,' on any file.')",mtime0
  return
endif
if (ifile1 eq -1) then begin
  print,format="('>>> Could not find end mtime ',3i4,' on any file.')",mtime1
  return
endif
;
; First file loop: get coord arrays and total number of times:
;
ntimes = 0
dtime1 = 0.
for ifile=ifile0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  finfo = ncdf_inquire(ncids[ifile])
  ndims = finfo.ndims
;
; Get spatial coord array dims from first file:
  if (ifile eq ifile0) then begin
    info1 = finfo
    ndims1 = info1.ndims
    dims1 = intarr(ndims1)
    for i=0,ndims1-1 do begin
      ncdf_diminq,ncids[ifile],i,name,len
      if (name eq 'lon') then nlon1=len
      if (name eq 'lat') then nlat1=len
      if (name eq 'lev') then nlev1=len
    endfor
  endif
;
; Coord array dims must be same as first file:
  for i=0,ndims-1 do begin
    ncdf_diminq,ncids[ifile],i,name,len
    for i=0,ndims-1 do begin
      ncdf_diminq,ncids[ifile],i,name,len
      if (name eq 'lon' and len ne nlon1) then begin
        print,'>>> lon dim (file 1)=',nlon1,' but lon dim (file ',ifile,')=',len
        return
      endif
      if (name eq 'lat' and len ne nlat1) then begin
        print,'>>> lat dim (file 1)=',nlat1,' but lat dim (file ',ifile,')=',len
        return
      endif
      if (name eq 'lev' and len ne nlev1) then begin
        print,'>>> lev dim (file 1)=',nlev1,' but lev dim (file ',ifile,')=',len
        return
      endif
    endfor
  endfor 
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    case vinfo.name of
      "mtime": ncdf_varget,ncids[ifile],ivar,mtime
      "lon"  : ncdf_varget,ncids[ifile],ivar,lon
      "lat"  : ncdf_varget,ncids[ifile],ivar,lat
      "lev"  : begin
        ncdf_varget,ncids[ifile],ivar,lev
        ncdf_attget,ncids[ifile],ivar,'units',units
        lev_units = string(units)
      end
      "time" : begin
        ncdf_varget,ncids[ifile],ivar,time
        ntime = n_elements(time)
        if (ifile eq ifile0) then ntime = ntime-ihist0
        if (ifile eq ifile1) then ntime = ihist1+1
        if (ifile0 eq ifile1 and ifile eq ifile0) then ntime = ihist1-ihist0+1
      end
      else:
    endcase
  endfor ; for ivar=0,nvars-1
  nlon = n_elements(lon)
  nlat = n_elements(lat)
  nlev = n_elements(lev)
  ntimes = ntimes+ntime
;
; Check for constant delta-time:
;
; if (ntime gt 0 and dtime1 eq 0.) then dtime1 = time[1]-time[0]
; for i=1,n_elements(time)-1 do begin
;   dtime = time[i]-time[i-1]
;   if (dtime ne dtime1) then begin
;     print,'>>> Delta-time has changed: dtime1=',dtime1,' current i=',i,$
;       ' time[i]=',time[i],' time[i-1]=',time[i-1],' dtime=',dtime
;     print,'Current file: ',files[ifile],' (ifile=',ifile,')'
;     return
;   endif
; endfor
  ncdf_attget,ncids[ifile],'missing_value',missing,/global
  ncdf_close,ncids[ifile]
  if (ifile eq ifile1) then break
endfor ; ifile=ifile0,nfiles-1
print,format="('nlon=',i4,' nlat=',i4,' nlev=',i4,' total ntimes=',i4)",$
  nlon,nlat,nlev,ntimes
;
; Use nearest grid-point to each requested lat,lon location:
;
ilocs = intarr(2,nlocs) ; indices to nearest grid point for each location
gridlocs = fltarr(2,nlocs) ; nearest grid-point lat,lon for each location
for iloc=0,nlocs-1 do begin
  ilat = ixfind(lat,loc_lats[iloc],lat[1]-lat[0])
  ilon = ixfind(lon,loc_lons[iloc],lon[1]-lon[0])
  ilocs[0,iloc] = ilat
  ilocs[1,iloc] = ilon
  gridlocs(0,iloc) = [lat[ilat]] ; actual lat grid point
  gridlocs(1,iloc) = [lon[ilon]] ; actual lon grid point
endfor
;
; Set up zp array for y-axis:
;
dzp = lev[1]-lev[0] ; assume uniform vertical grid
if (zpbot eq zptop) then begin ; default to full zp grid
  zplev = lev
  nzp = nlev
  k0 = 0
  k1 = nlev-1
endif else begin
  if (zpbot lt lev[0]) then begin ; default to bottom of grid
    k0 = 0
  endif else begin
    k0 = ixfind(lev,zpbot,dzp)
    if (k0 eq -1) then begin
      print,'>>> Could not find grid point for zpbot=',zpbot
      print,'levels = ',lev
      return
    endif
  endelse
  if (zptop gt lev[nlev-1]) then begin ; default to top of grid
    k1 = nlev-1
  endif else begin
    k1 = ixfind(lev,zptop,dzp)
    if (k1 eq -1) then begin
      print,'>>> Could not find grid point for zptop=',zptop
      print,'levels = ',lev
      return
    endif
  endelse
  nzp = k1-k0+1
  zplev = fltarr(nzp)
  zplev[0:nzp-1] = lev[k0:k1]
endelse
print,format="('nzp=',i3,' k0,k1=',i3,',',i3,' zplev=',f9.3,' to ',f9.3)",$
  nzp,k0,k1,zplev[0],zplev[nzp-1]
;
; Main field array:
;
f = dblarr(nlocs,nzp,ntimes,nfields)
mtimes = intarr(3,ntimes)
years = intarr(ntimes)    ; integer years
decyears = fltarr(ntimes) ; decimal years
itime0 = 0
nyrbound = 0
iyrbound = intarr(20) ; assume max of 20 years
;
; File loop to read vars:
;
for ifile=ifile0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  finfo = ncdf_inquire(ncids[ifile])
;
; Get number of times (histories) on current file and update: 
;   ntime: Number of times in current file.
;   time0,time1: Beginning and ending time indices into 
;     total-time arrays (0->ntimes-1) for current file.
;
  ndims = finfo.ndims
  for idim=0,ndims-1 do begin
    ncdf_diminq,ncids[ifile],idim,dimname,dimsize
    if (dimname eq 'time') then begin
      ntime = dimsize 
;
; Number of histories to read from this file:
      if (ifile eq ifile0) then ntime = ntime-ihist0
      if (ifile eq ifile1) then ntime = ihist1+1
      if (ifile0 eq ifile1 and ifile eq ifile0) then ntime = ihist1-ihist0+1
      itime1 = itime0+ntime-1 
    endif
  endfor
;
; First,last histories to read from this file:
;
  ihistbeg = 0 & ihistend = ntime-1
  if (ifile eq ifile0) then begin
    ihistbeg = ihist0 & ihistend = ihist0+ntime-1
  endif else if (ifile eq ifile1) then begin
    ihistbeg = 0 & ihistend = ihist1
  endif
; print,format="(/,'itime0,1=',2i4,' ihistbeg,end=',2i4,' ntime=',i4)",$
;   itime0,itime1,ihistbeg,ihistend,ntime
;
; Get years:
;
  id = ncdf_varid(ncids[ifile],'year')
  if (id eq -1) then print,'>>> WARNING: could not find year variable'
  ncdf_varget,ncids[ifile],id,var
; print,' ' & print,'itime0=',itime0,' itime1=',itime1,' ntime=',ntime
  years[itime0:itime1] = fix(var[ihistbeg:ihistend])
;
; Get model times:
;
  id = ncdf_varid(ncids[ifile],'mtime')
  if (id eq -1) then print,'>>> WARNING: could not find mtime variable'
  ncdf_varget,ncids[ifile],id,var
  mtimes[*,itime0:itime1] = var[*,ihistbeg:ihistend]
;
; Check for year-boundary:
;
  for i=itime0,itime1 do begin
    decyears[i] = float(years[i])+float(mtimes[0,i])/365.
    if (decyears[i] eq float(fix(decyears[i]))) then begin
      print,'Found year-boundary at i=',i,' decyears=',decyears[i]
      iyrbound[nyrbound] = i
      nyrbound = nyrbound+1
    endif
  endfor
;
; Print time info for current file:
;
  print,format="(/,'File ',a,' (',i4,' of ',i4,'):')",$
    file_basename(files[ifile]),ifile-ifile0+1,ifile1-ifile0+1
  print,format=$
    "('Reading ',i3,' model times: ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3)",$
    ntime,mtimes[*,itime0],mtimes[*,itime1]
;
; Read requested 4d fields:
;
  found = intarr(nfields)
  found[*] = 0
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    for ifld=0,nfields-1 do begin
      if (vinfo.name eq fields[ifld]) then begin
        found[ifld] = 1
        if (vinfo.ndims ne 4) then begin ; should also verify dim lens
          print,'>>> Field ',fields[ifld],' is not 4-d: skipping..'
          f[*,*,itime0:itime1,ifld] = 0.
          break
        endif
        ncdf_varget,ncids[ifile],ivar,var ; (lon,lat,lev,time)
        print,format="('Field ',a,': 4d min,max=',2e12.4)",fields[ifld],$
          min(var[*,*,k0:k1,ihistbeg:ihistend]),max(var[*,*,k0:k1,ihistbeg:ihistend])
        for iloc=0,nlocs-1 do begin
;
; f = fltarr(nlocs,nlev,ntimes,nfields)
          f[iloc,*,itime0:itime1,ifld] = $
            var[ilocs[1,iloc],ilocs[0,iloc],k0:k1,ihistbeg:ihistend]
          if (k1 eq nlev-1) then $
            f[iloc,nzp-1,itime0:itime1,ifld] = f[iloc,nzp-2,itime0:itime1,ifld]
        endfor
;
; Check for units attribute of current variable (may not exist for all fields):
        f_units[ifld] = ""
        for iatt=0,vinfo.natts-1 do begin
          attname = ncdf_attname(ncids[ifile],ivar,iatt)
          if (attname eq "units") then begin
            ncdf_attget,ncids[ifile],ivar,attname,units
            f_units[ifld]=string(units)
          endif
        endfor
      endif ; is a requested field
    endfor ; requested field loop
  endfor ; ivar=0,nvars-1
  for ivar=0,nfields-1 do begin
    if (found[ivar] eq 0) then print,format=$
      "('(Field ',a,' was not found)')",fields[ivar]
  endfor
  itime0 = itime0+ntime
  ncdf_close,ncids[ifile]
  if (ifile eq ifile1) then break
endfor ; ifile=ifile0,nfiles-1
if (ntimes le 1) then begin
  print,'>>> ntimes=',ntimes,': must have more than 1 time.'  
  return
endif
;
; Insure consistent delta-time (minutes) between all histories.
;
; If a deltamin between 2 histories is zero (duplicate history), 
; then the arrays are shifted to overwrite the 1st history of the 
; duplicate pair, and ntimes is decremented.  Duplicate secondary 
; histories can occur between the end of one file and the beginning 
; of the next when the model is restarted.  If deltamin is inconsistent 
; w/ deltamin0 and deltamin is non-zero, then this is a fatal
; (the history write frequency was apparently changed mid-run).
;
deltamin0 = mtime_to_mins(mtimes[*,1])-mtime_to_mins(mtimes[*,0])
dtimeloop:
for i=1,ntimes-1 do begin
  deltamin = mtime_to_mins(mtimes[*,i])-mtime_to_mins(mtimes[*,i-1])
  if (deltamin ne deltamin0) then begin
    print,format="(/,'Warning: deltamin is inconsistent: deltamin0=',i5,' deltamin=',i5)",$
      deltamin0,deltamin
    print,format="('      i=',i5,' mtimes[i]=',3i4,' mtimes[i-1]=',3i4)",$
      i,mtimes[*,i],mtimes[*,i-1]
    if (deltamin ne 0) then return
    print,format="('      (Duplicate histories will be corrected by shifting the arrays)')"
    mtimes[*,i:ntimes-1] = shift(mtimes[*,i:ntimes-1],0,-1)
    years[i:ntimes-1] = shift(years[i:ntimes-1],-1)
    decyears[i:ntimes-1] = shift(decyears[i:ntimes-1],-1)
    for ii=0,nfields-1 do begin
      f[*,*,i:ntimes-1,ii] = shift(f[*,*,i:ntimes-1,ii],0,0,-1)
    endfor
; This results in an error in the number of arguments to shift function:
;   f[*,*,i:ntimes-1,*] = shift(f[*,*,i:ntimes-1,*],0,0,-1,0)
    ntimes = ntimes-1
    goto,dtimeloop ; check for another one
  endif
endfor
;
print,format=$
  "(/,'Total ntimes=',i5,': model times ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3,' by ',i6,' minutes')",$
  ntimes,mtimes[*,0],mtimes[*,ntimes-1],deltamin
;
; Set plot position:
;
xll = .10 & yll = .33
xur = .95 & yur = yll+.7*(xur-xll)
axespos  = [xll,yll,xur,yur]
xmid = .5*(xur+xll)
;
; Graphics device is X or ps:
;
if (keyword_set(ps)) then begin
  set_plot,'PS'
  device,/color,/landscape                             ; default width=9.5" (prints on 8.5x11)
; device,/color,/landscape,/inches,xsize=11.,ysize=8.  ; 11" wide
; device,/color,/landscape,/inches,xsize=14.,ysize=8.  ; fits 1 yr w/ scrolling
; device,/color,/landscape,/inches,xsize=25.,ysize=8.  ; best for > 2-3 years
; device,/color,/landscape,/inches,xsize=35.,ysize=8.
endif else begin
  set_plot,'X'
  device,decomposed=0,retain=2 ; for TrueColor display
endelse

 nlevels = 12 ; number of contour levels
 loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red
;
; For horizontal color bar at bottom:
;
barw = xur-xll-.05*(xur-xll) & barh = .03
xbar = xmid-.5*barw
ybar = yll-barh-.12
baroff = .035
;
;chsize = 1.2
 chsize = 1.0
;
; Field plot loop:
;
iframe = 0
for ifld=0,nfields-1 do begin
  if (found[ifld] eq 0) then begin
    print,'Skipping field ',fields[ifld],' because it was not found.'
    continue ; field was not found
  endif
  for iloc=0,nlocs-1 do begin
    utv = f[iloc,*,0:ntimes-1,ifld] ; [nlev,ntimes]
    utv = transpose(utv)   ; [ntimes,nlev]
    funits = f_units[ifld]
    if (f_units[ifld] eq 'CM') then begin
      utv = utv/1.e5          ; cm to km (e.g., Z)
      funits = 'KM'
    endif
    if (f_units[ifld] eq 'CM/S') then begin
      utv = utv/100.
      funits = 'M/S'
    endif
;
; Contour levels:
;
    fmin = min(utv(where(utv ne missing)))
    fmax = max(utv(where(utv ne missing)))
    step = (fmax-fmin)/nlevels
    if (step le 0.) then begin
      print,format=$
        "('>>> Contour interval=',e12.4,' -- skipping field ',a,' at lat,lon=',f6.2,',',f7.2)",$
        step,fields[ifld],loc_lats[iloc],loc_lons[iloc]
      print,format="('    fmin,fmax=',2e12.4)",fmin,fmax
      continue
    endif
    clevels = indgen(nlevels)*step+fmin
    c_colors = indgen(nlevels)+3
;
; Main title and y-axis title:
;
    title = fields[ifld]+'('+funits+'):  Lat,Lon = '+ $
      string(format="(f6.2)",gridlocs[0,iloc])+', '+$
      string(format="(f7.2)",gridlocs[1,iloc])
    ytitle = "Log Pressure ("+lev_units+")"
;
; Use IDL procedures timegen, julday, and label_date 
; (setting nxticks does not appear to work w/ label_date)
;
    time = timeaxis(years[0:ntimes-1],mtimes[*,0:ntimes-1],deltamin,xtitle,$
      nxticks)
;
; First contour call: save xtick values with xtick_get, but do not draw xaxes 
; or data yet.
;
    contour,utv,time,zplev,position=axespos,/nodata,$
      xstyle=5,ystyle=5,xtick_get=xtickv

;     xstyle=5,ystyle=5,xticks=nxticks,xtick_get=xtickv
;
; This uses the time coord var calculated in timeaxis function 
; to make major x-axis tick labels month/day/year.
    date_label = label_date(date_format="%N/%D/%Y!C")
;
; Set number of minor tick marks:
;
    for i=0,1 do begin
      caldat,xtickv[i],calmon,calday,calyear,calhr,calmin,calsec
      if (i gt 0) then dticks = xtickv[i]-xtickv[i-1]
    endfor
    if (dticks eq 1.0) then begin ; 1 days between major ticks 
      xminor = 12
    endif else if (dticks eq 0.5) then begin ; 12 hours between major ticks
      xminor = 12
    endif else if (dticks eq 1./24.) then begin ; 1 hour between major ticks
      xminor = 6
    endif else xminor = 0
;
; Second contour call with data and axes:
;
; Use IDL pro label_date from above:
    contour,utv,time,zplev,position=axespos,/noerase,$
      /fill,ytitle=ytitle,title=title,xminor=xminor,$
      charsize=chsize,min_value=fmin,$
      xstyle=1,ystyle=1,xtickformat='LABEL_DATE',$
      c_colors=c_colors,levels=clevels,background=255,$
      color=1,max_value=fmax
;
; Add local time to xtick labels below date label:
; (note 0.5 days are added to the julian day xtick value since
;  julian days start at 12 noon (s.a., timeaxis function above)
;
    nxtick = n_elements(xtickv)
    for i=0,nxtick-1 do begin
      caldat,xtickv[i]+0.5,calmon,calday,calyear,calhr,calmin,calsec
      uthr = float(calhr)+float(calmin)/60.
      slt = calcslt(uthr,loc_lons[iloc])
      lthr = fix(slt)
      ltmin = fix((slt-float(lthr))*60.)
;
; Get day of year:
; function monthday_to_dayofyear, imonth,iday
      dayofyear = monthday_to_dayofyear(calmon,calday)
      dayofyear = dayofyear-1 ; because the julian day starts at noon
;
      sltlab = string(format="('day ',i3.3,' lt ',i2.2,':',i2.2)",$
        dayofyear,lthr,ltmin)
      xnorm = convert_coord(xtickv[i],0.,/to_normal)
;     xyouts,xnorm[0,0],yll-.055,sltlab,/norm,align=0.5,charsize=chsize-0.2
      xyouts,xnorm[0,0],yll-.050,sltlab,/norm,align=0.5,charsize=chsize-0.2
    endfor
;
; X-axis label
    xyouts,xmid,yll-.09,xtitle,/norm,align=0.5,charsize=chsize
;
; Place vertical dashed line at any year-boundaries:
;
    if (nyrbound gt 0) then begin
      for i=0,nyrbound-1 do begin
        oplot,[iyrbound[i],iyrbound[i]],[zplev[0],zplev[nzp-1]],linestyle=2 ; dashed line
      endfor
    endif
;
; Min,max label:
;
    cint = clevels[1]-clevels[0] ; assume constant contour interval
    minmaxlabel = string(format="('Min,Max = ',g10.4,', ',g10.4,' Interval = ',g10.4)",$
      fmin,fmax,cint)
    xyouts,xmid,yll-.24,minmaxlabel,/norm,align=0.5,charsize=chsize
;
; Files label:
;
    filelabel = "Files " + file0 + " to " + file1
    xyouts,xmid,yll-.28,filelabel,/norm,align=0.5,charsize=chsize
;
; Color bar:
;
    form = "(g10.4)"
;   form = "(f9.3)"
    cblabels = strarr(nlevels+1)
    cblabels[*] = ' '
    cblabels[0] = strcompress(string(format=form,fmin),/remove_all)
    for i=2,nlevels-2,2 do begin ; label every other level on color bar
      cblabels[i] = strcompress(string(format=form,clevels[i]),/remove_all)
    endfor
    cblabels[nlevels] = strcompress(string(format=form,fmax),$
      /remove_all)
    clrbar,barw,barh,xbar,ybar,cblabels,botcolor=3,topcolor=nlevels+3-1,$
      charsize=chsize,laboffset=baroff,labcolor=1
    if (not keyword_set(ps)) then begin
      cursor,xcurs,ycurs,4
    endif
    iframe = iframe+1
    print,format="('Frame ',i3,' Field ',a,' lat,lon=',f6.2,', ',f7.2)",$
      iframe,fields[ifld],loc_lats[iloc],loc_lons[iloc]
  endfor ; iloc
endfor ; ifld
if (keyword_set(ps)) then begin
  device,/close
endif
end