;
; User input (called from main pro utproc below):
;
pro input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,zpbot,zptop
;
; Directories and filename roots of nc history files to read.
; (dirs and fileroots must have same number of elements)
;
 dirs = ["/oreo/d/foster/icnossen/sy1957ds"]
 fileroots = ["ICNOSSEN.tiegcm1.8.1957.sy1957ds_"]
;
; Beginning and ending file numbers to read:
;
 file_beg = [44]
 file_end = [47]
;
; Format of file numbers (2 or 3 digits):
;
;form = "(i02)"
 form = "(i03)"
;
; 4d fields to read and plot:
;
 fields = ["TN","UN","VN","O2","O1","NO","N4S","NE","TE","TI","O2P",$
           "W","Z","BE3","BTOT","SIN(I)","INCL"]
;fields = ["TN","UN","VN","NE","Z"]
;fields = ["TN","UN","VN","W","O1","O2","NO","N4S","NE","TE","TI","O2P",$
;          "POTEN","UTNDCY","VTNDCY","DIFKK","Z"]
;fields = ["TN","UN","VN","W","O1","O2","NO","N4S","NE","TE","TI","O2P","Z"]
;fields = ["TN","UN","VN","OMEGA","O1","O2","NO","N4S","NE","TE","TI","O2P","Z"]
;
; 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.]
 loc_lats = [  -80. , -70., -60., -50.,  50.,  60.,  70.,  80.]
;loc_lons = [-180.]
;loc_lats = [-62.5]
;
; 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 = -7.
;zptop = -1.
;zpbot = -7.
 zptop = 5.
;zpbot = 0. & zptop = 0.
;zpbot = -999. & zptop = 999.
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 mtime_to_mins,mtime
return,mtime(0)*24*60+mtime(1)*60+mtime(2)
end
;-----------------------------------------------------------------------
function timeaxis,years,mtimes,deltamin,xtitle
dayofyear_to_monthday,mtimes(0,0),imonth0,iday0
julian0 = julday(imonth0,iday0,years[0],mtimes(1,0),mtimes(2,0),0)
ntimes = n_elements(years)
time = timegen(ntimes,start=julian0,step_size=deltamin,units="Minutes")
;
; Daily:
;
if (deltamin eq 1440) then begin          ; daily
; date_label = label_date(date_format="%N/%D/%Y")
  date_label = label_date(date_format="%N/%D/%Z")
  if (years[ntimes-1]-years[0] le 0) then begin
    xtitle = string(format=$
      "('Time: Days ',i3,' to ',i3,' of ',i4,' (',i4,' points)')",$
      mtimes[0,0],mtimes[0,ntimes-1],years[0],ntimes) 
  endif else begin
    xtitle = string(format=$
      "('Time: Days ',i3,' of ',i4,' to day ',i3,' of ',i4,' (',i4,' points)')",$
      ntimes,mtimes[*,0],years[0],mtimes[*,ntimes-1],years[ntimes-1],ntimes)
  endelse
;
; Hourly:
;
endif else if (deltamin eq 60) then begin ; hourly
; date_label = label_date(date_format="%N/%D/%Z!C(hour%H)")
  date_label = label_date(date_format="%N/%D/%Z")
  if (years[ntimes-1]-years[0] le 0) then begin
    xtitle = string(format=$
      "('Time: Hourly from ddd:hh ',i3.3,':',i2,' to ',i3.3,':',i2,' of ',i4,' (',i4,' points)')",$
      mtimes[0,0],mtimes[1,0],mtimes[0,ntimes-1],mtimes[1,ntimes-1],years[0],ntimes)
  endif else begin
    xtitle = string(format=$
      "('Time: Hourly from yyyy:ddd:hh ',i4,':',i3.3,':',i2,' to ',i4,':',i3.3,':',i2,' (',i4,' points)')",$
      years[0],mtimes[0,0],mtimes[1,0], years[ntimes-1],mtimes[0,ntimes-1],mtimes[1,ntimes-1],ntimes)
  endelse
;
; Generic (not daily or hourly)
;
endif else begin                          ; generic
; date_label = label_date(date_format="%N/%D/%Y!C(%H:%I)")
  date_label = label_date(date_format="%N/%D/%Y")
  xtitle = string(format=$
    "('TIME (',i5,' points from year:day:hour:min ',i4,':',i3.3,':',i2.2,':',i2.2,' to ',i4,':',i3.3,':',i2.2,':',i2.2,' by ',i4,' minutes)')",$
    ntimes,years[0],mtimes[*,0],years[ntimes-1],mtimes[*,ntimes-1],deltamin)
endelse
time = time+0.5 ; julian days start at hour 12, tgcm days start at hr 0
return,time
end
;-----------------------------------------------------------------------
; Main:
;
pro utproc,ps=ps
input,dirs,fileroots,file_beg,file_end,form,fields,loc_lons,loc_lats,zpbot,zptop
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
ntimes = 0
;
; First file loop: get coord arrays and total number of times:
;
for ifile=0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  finfo = ncdf_inquire(ncids[ifile])
  ndims = finfo.ndims
  if (ifile eq 0) then begin
    info1 = finfo
    ndims1 = info1.ndims
    dims1 = intarr(ndims1)
    for i=0,ndims1-1 do begin
      ncdf_diminq,ncids[ifile],i,name,len
      dims1[i] = len
    endfor
  endif
  if (ndims ne ndims1) then $
    print,format="('>>> ndims=',i3,' ne ndims1=',i3)",ndims,ndims1
  dims = intarr(ndims)
  for i=0,ndims-1 do begin
    ncdf_diminq,ncids[ifile],i,name,len
    dims[i] = len 
    if (i ne finfo.recdim and dims[i] ne dims1[i]) then begin
      print,format="('>>> i=',i3,' dims1[i]=',i4,' ne dims[i]=',i4)",$
        i,dims1[i],dims[i]
      return
    endif
  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
      else:
    endcase
  endfor
  ntime = n_elements(mtime)/3
  nlon = n_elements(lon)
  nlat = n_elements(lat)
  nlev = n_elements(lev)
  ntimes = ntimes+ntime
  ncdf_attget,ncids[ifile],'missing_value',missing,/global
  ncdf_close,ncids[ifile]
endfor ; ifile=0,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 = fltarr(nlocs,nzp,ntimes,nfields)
mtimes = intarr(3,ntimes)
ixindex = lonarr(ntimes) ; long-integer array
for i=0,ntimes-1 do ixindex[i]=i
years = intarr(ntimes)    ; integer years
decyears = fltarr(ntimes) ; decimal years
itime0 = 0
nyrbound = 0
iyrbound = intarr(20) ; assume max of 20 years
for ifile=0,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 
      itime1 = itime0+ntime-1 
    endif
  endfor
;
; 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
  years[itime0:itime1] = fix(var)
;
; 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
;
; 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 ',i4,' of ',i4,': ',a)",ifile+1,nfiles,files[ifile]
  print,format=$
    "('Reading model times ',i3,',',i3,',',i3,' to ',i3,',',i3,',',i3)",$
    mtimes[*,itime0],mtimes[*,itime1]
;
; Read requested 4d fields:
;
  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
        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),max(var)
        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,*]
          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
  itime0 = itime0+ntime
  ncdf_close,ncids[ifile]
endfor ; ifile=0,nfiles-1
;
; Determine delta time between histories and insure it is constant:
;
 if (ntimes le 1) then begin
   print,'>>> ntimes=',ntimes,': must have more than 1 time.'  
   return
 endif
 deltamin0 = mtime_to_mins(mtimes[*,1])-mtime_to_mins(mtimes[*,0])
 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,'>>> deltamin is inconsistent: deltamin0=',deltamin0,$
       ' i=',i,' deltamin=',deltamin
   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 size
; device,/color,/landscape,/inches,xsize=14.,ysize=8.  ; fits up to 1 yr w/o 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
  for iloc=0,nlocs-1 do begin
    utv = f[iloc,*,*,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
    fmin = min(utv)
    fmax = max(utv)
;
; Contour levels:
;
    step = (fmax-fmin)/nlevels
    if (step le 0.) then begin
      print,'>>> bad contour interval step=',step,' -- skipping this plot'
      break
    endif
    clevels = indgen(nlevels)*step+fmin
;   print,format=$
;     "('number of contour levels=',i3,' fmin,max=',2e12.4,' step=',e12.4,' clevels=',/,(6e12.4))",$
;     nlevels,fmin,fmax,step,clevels 
    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 
;
    time = timeaxis(years,mtimes,deltamin,xtitle)
;   for i=0,n_elements(time)-1 do begin
;     caldat,time[i],calmon,calday,calyear,calhr,calmin,calsec
;     print,format=$
;       "('i=',i4,' julian time[i]=',f10.2,' mon=',f6.2,' day=',f6.2,' year=',f8.2,' hr,min,sec=',3f9.2)",$
;       i,time[i],calmon,calday,calyear,calhr,calmin,calsec
;   endfor
;
; 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

;   for i=0,n_elements(xtickv)-1 do begin
;     caldat,xtickv[i],calmon,calday,calyear,calhr,calmin,calsec
;     print,format=$
;       "('major xticks: i=',i4,' julian xtickv[i]=',f10.2,' mon=',f6.2,' day=',f6.2,' year=',f8.2,' hr,min,sec=',3f9.2)",$
;       i,xtickv[i],calmon,calday,calyear,calhr,calmin,calsec
;   endfor
    xminor = 12
;
; 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,max_value=fmax,$
      xstyle=1,ystyle=1,xtickformat='LABEL_DATE',$
      c_colors=c_colors,levels=clevels,background=255,color=1
;
; X-axis label
    xyouts,xmid,yll-.08,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