;-----------------------------------------------------------------------
; 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 utproc,ps=ps
;
; Directories containing nc history files to read:
;
;dir = "/hao/tgcm/timegcm1.2_y92ur/py92ur"     ; cisl
;dirs = ["/aim/d/tgcm/timegcm1.2_y92ur/py92ur"]   ; hao
;
 dirs = ["/aim/d/tgcm/timegcm1.2_y92ur/py92ur",$ 
         "/aim/d/tgcm/timegcm1.2_y93ur/py93ur",$
         "/aim/d/tgcm/timegcm1.2_y94ur/py94ur"]
;
 fileroots = ["ROBLE.timegcm1.2.py92ur",$   ; root of file names in each dir
              "ROBLE.timegcm1.2.py93ur",$
              "ROBLE.timegcm1.2.py94ur"]
;
;file_beg = [1,1]    ; beginning file number for each dir
;file_end = [37,37]  ; ending file number for each dir
;
;file_beg = [35,1] ; 55 times
;file_end = [37,3]

;file_beg = [20,1] ; 275 times
;file_end = [37,10]

 file_beg = [1,1,1]
 file_end = [37,37,37]
;
 form = "(i02)"
;fields = ['O2']
;fields = ['TN']
;fields = ['TN','UN','VN','O2','Z']
 fields = ["TN","UN","VN","O2","O1","O3","H2O","CO2","NE","Z"]
 loc_lons = [-180.,-180.,-180.,-180.,-180.]
 loc_lats = [ -60., -30.,   0.,  30.,  60.]
;loc_lons = [  0.]
;loc_lats = [-60.]
;
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)
    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,'nfiles=',nfiles,' nfields=',nfields,' nlocs=',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"  : ncdf_varget,ncids[ifile],ivar,lev
      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,'nlon=',nlon,' nlat=',nlat,' nlev=',nlev,' total ntimes=',ntimes
;
; Main field array:
;
f = fltarr(nlocs,nlev,ntimes,nfields)
mtimes = intarr(3,ntimes)
days = intarr(ntimes)
ixindex = lonarr(ntimes) ; long-integer array
for i=0,ntimes-1 do ixindex[i]=i
iyyyyddd = lonarr(ntimes) ; long-integer array
years = fltarr(ntimes)
itime = 0
for ifile=0,nfiles-1 do begin
  ncids[ifile] = ncdf_open(files[ifile],/nowrite)
  print,format="('File ',i3,': ',a,' ncid=',i3)",ifile,files[ifile],ncids[ifile]
  finfo = ncdf_inquire(ncids[ifile])
;
; Get years:
;
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    if (vinfo.name eq "year") then begin
      ncdf_varget,ncids[ifile],ivar,var
      ntime = n_elements(var)
      itime1 = itime+ntime-1
      years[itime:itime1] = var
      break 
    endif
  endfor ; ivar for iyear
;
; Calculate yyddd: 
;
  for ivar = 0,finfo.nvars-1 do begin
    vinfo = ncdf_varinq(ncids[ifile],ivar)
    if (vinfo.name eq "mtime") then begin
      ncdf_varget,ncids[ifile],ivar,var
      ntime = n_elements(var)/3
      itime1 = itime+ntime-1
      mtimes[*,itime:itime1] = var
      for i=itime,itime1 do $
        iyyyyddd[i] = years[i]*1000.+mtimes[0,i]
      break
    endif ; mtime var
  endfor ; ivar for mtime
;
; Read fields to plot:
;
  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
        ncdf_varget,ncids[ifile],ivar,var
        if (vinfo.ndims ne 4) then begin ; should also verify dim lens
          print,'>>> Field ',fields[ifld],' is not 4-d: skipping..'
          f[*,*,itime:itime1,ifld] = 0.
          break
        endif
        for iloc=0,nlocs-1 do begin
          ilon = ixfind(lon,loc_lons[iloc],lon[1]-lon[0])
          ilat = ixfind(lat,loc_lats[iloc],lat[1]-lat[0])
;
; f = fltarr(nlocs,nlev,ntimes,nfields)
          f[iloc,*,itime:itime1,ifld] = var[ilon,ilat,*,*]
          f[iloc,nlev-1,itime:itime1,ifld] = f[iloc,nlev-2,itime:itime1,ifld]
        endfor
        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
            units = string(units)
            f_units[ifld]=units
          endif
        endfor
        print,format=$
          "('Read field ',a,' times ',i4,' to ',i4,' of ',i4)",$
          fields[ifld],itime,itime1,ntimes
      endif
    endfor
  endfor ; ivar=0,nvars-1
  itime = itime+ntime
  ncdf_close,ncids[ifile]
endfor ; ifile=0,nfiles-1
;
; 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
; device,/color,/inches,xsize=10.,ysize=5.
; device,/color,/landscape,/inches,xsize=15.,ysize=8.
; device,/color,/landscape,/inches,xsize=25.,ysize=8.
  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
;
; Set up x-axis:
;
  nxticks = 8
  xtickspace = ntimes/(nxticks-1)
  nxminor = xtickspace
  for i=10,100,5 do begin
    if (xtickspace mod i eq 0) then nxminor = i
  endfor
  print,'nxticks=',nxticks,' xtickspace=',xtickspace,' nxminor=',nxminor
  xticknames = strarr(nxticks)
  xtickv = lonarr(nxticks)
  itick = 0
  for i=0,ntimes-1 do begin
    if (i mod xtickspace eq 0) then begin
      xticknames[itick] = string(format="(i7)",iyyyyddd[i])
      xtickv[itick] = ixindex[i]
      itick = itick+1
    endif
  endfor
;
; 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.1
;
; 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]
    if (f_units[ifld] eq 'CM') then begin
      utv = utv/1.e5
      f_units[ifld] = 'KM'
    endif
    if (f_units[ifld] eq 'CM/S') then begin
      utv = utv/1.e3
      f_units[ifld] = 'M/S'
    endif
    fmin = min(utv)
    fmax = max(utv)
;
; Contour levels:
;
    step = (fmax-fmin)/nlevels
    clevels = indgen(nlevels)*step+fmin
    c_colors = indgen(nlevels)+3
;
; Title labels:
;
    title = fields[ifld]+'('+f_units[ifld]+'):  Lon,Lat = '+ $
      string(format="(f7.2)",loc_lons[iloc])+', '+ $
      string(format="(f6.2)",loc_lats[iloc])
    xtitle = string(format="('YearDay (',i4,' days from ',i7,' to ',i7,')')",$
      ntimes,iyyyyddd[0],iyyyyddd[ntimes-1])
    ytitle = "Zp"
;
; Contour:
;
    contour,utv,ixindex,lev,xstyle=1,ystyle=1,position=axespos,$
      /fill,xtitle=xtitle,ytitle=ytitle,title=title,$
      charsize=chsize,min_value=fmin,max_value=fmax,$
      xticks=nxticks,xtickv=xtickv,xtickname=xticknames,xminor=nxminor,$
      c_colors=c_colors,levels=clevels,background=255,color=1
;
; Min,max label:
;
    minmaxlabel = string(format="('Min,Max = ',g9.3,', ',g9.3,' (',a,')')",$
      fmin,fmax,f_units[ifld])
    xyouts,xmid,yll-.24,minmaxlabel,/norm,align=0.5,charsize=chsize
;
; Color bar:
;
    cblabels = strarr(nlevels+1)
    cblabels[*] = ' '
    cblabels[0] = strcompress(string(format="(g8.2)",fmin),/remove_all)
    for i=3,nlevels-3,3 do begin ; label every 3rd level on color bar
      cblabels[i] = strcompress(string(format="(g8.2)",clevels[i]),/remove_all)
    endfor
    cblabels[nlevels] = strcompress(string(format="(g8.2)",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