;-----------------------------------------------------------------------
function mkclevels,fmin,fmax,cmin,cmax,cint
;
; Set up contour min,max,interval. Args fmin,fmax are input,
; args cmin,cmax,cint are intent inout.
;
;print,'Enter mkclevels: fmin,max=',fmin,fmax,' cmin,max=',cmin,cmax,$
;  ' cint=',cint
;
; Check for valid fmin,fmax (e.g., all data missing)
;
if fmin ge fmax then begin
  print,'>>> WARNING mkclevels: fmin ge fmax: fmin=',fmin,' fmax=',fmax
  return,[0.]
endif
;
; Default number of levels:
;
;nlevdef = 12
nlevdef = 10
;
; If cmin >= cmax and cint <= 0, do not use cmin,cmax, or cint:
;
if cmin ge cmax and cint le 0. then begin
  nlevels = nlevdef
  levels = fltarr(nlevels)
  conint = (fmax-fmin) / (nlevels-1)
  for i=0,nlevels-1 do levels(i) = fmin+i*conint
  if cmin eq cmax then begin
    cmin=fmin & cmax=fmax
  endif
  cint = conint
; print,'mkclevels: cmin,max,int=',cmin,cmax,cint,' levels=' & print,levels
  return,levels
endif
;
; If cmin < cmax but cint <= 0, then use cmin and cmax, but not cint:
;
if cmin lt cmax and cint le 0. then begin
  nlevels = nlevdef
  levels = fltarr(nlevels)
  cint = (cmax-cmin) / (nlevels-1)
  for i=0,nlevels-1 do levels(i) = cmin+i*cint
; print,'mkclevels: cmin,max,int=',cmin,cmax,cint,' levels=' & print,levels
  return,levels
endif
;
; If cmin >= cmax and cint > 0, then use cint but not cmin or cmax:
;
if cmin ge cmax and cint gt 0. then begin
  nlevels = fix((fmax-fmin)/cint+.5) ; +.5 is to get nearest int
  if nlevels gt 59 then begin
    cintold = cint
    while nlevels gt 59 do begin
      cint = cint+.01*cint
      nlevels = fix((fmax-fmin)/cint+.5) ; +.5 is to get nearest int
    end
    print,'WARNING: had to increase cint from ',cintold,' to ',cint,$
      ' to avoid > 60 levels'
  endif
  if cmin+float(nlevels)*cint le cmax then nlevels = nlevels+1
  levels = fltarr(nlevels)
  for i=0,nlevels-1 do levels(i) = fmin+i*cint
  if cmin eq cmax then begin
    cmin = fmin & cmax = fmax
  endif
; print,'mkclevels: cmin,max,int=',cmin,cmax,cint,' levels=' & print,levels
  return,levels
endif
;
; If cmin < cmax and cint > 0, then use cmin, cmax, and cint:
;
if cmin lt cmax and cint gt 0. then begin
  nlevels = fix((cmax-cmin)/cint+.5) ; +.5 is to get nearest int
  if nlevels gt 59 then begin
    cintold = cint
    while nlevels gt 59 do begin
      cint = cint+.01*cint
      nlevels = fix((fmax-fmin)/cint+.5) ; +.5 is to get nearest int
    end
    print,'WARNING: had to increase cint from ',cintold,' to ',cint,$
      ' to avoid > 60 levels: nlevels now =',nlevels
  endif
  if cmin+float(nlevels)*cint le cmax then nlevels = nlevels+1
  levels = fltarr(nlevels)
  for i=0,nlevels-1 do levels(i) = cmin+i*cint
; print,'mkclevels: cmin,max,int=',cmin,cmax,cint,' levels=' & print,levels
  return,levels
endif
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 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
;
;-----------------------------------------------------------------------
; 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
;
;-----------------------------------------------------------------------
;
pro utcontour,stateptr,f,time,yy,glon,cmin,cmax,cint,xtitle,ytitle,title,$
  filelabel,missing,jpg=jpg,ps=ps,png=png

  state = *stateptr
;
; Set plot position:
;
  xll = .10 & yll = .33
  xur = .95 & yur = yll+.7*(xur-xll)
  axespos  = [xll,yll,xur,yur]
  xmid = .5*(xur+xll)
;
; 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
;
; Contour levels:
;
  indices = where(f ne missing)
  if (n_elements(indices) eq 1) then begin
    if (indices eq -1) then begin
      print,'>>> All values are missing (no plot made): ',title
      title = '' ; this is a flag to the calling procedure
      return
    endif
  endif
;
; fmin,max = field min,max where not missing:
; cmin,cmax = optionally user supplied (default = 0.,0.)
;
  fmin = min(f(where(f ne missing)))
  fmax = max(f(where(f ne missing)))
  if (cmin ge cmax) then begin ; user did not provide cmin,cmax
    cmin = fmin
    cmax = fmax
  endif
;
; Function mkclevels from tgcmproc_idl:
;
  clevels = mkclevels(fmin,fmax,cmin,cmax,cint)
  nlevels = n_elements(clevels)
; print,'utcontour: fmin,fmax=',fmin,',',fmax,' cmin,max,int=',$
;   cmin,',',cmax,',',cint,' nlevels=',nlevels,' clevels=',clevels

  loadct,33,ncolors=nlevels,bottom=3 ; Blue-Red

  c_colors = indgen(nlevels)+3
;
; contour_label = 0 ; options of labeling contours
; contour_label = 1 ; options of labeling contours
  contour_label = state.label_contours
;
; 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")
;
; Get xtickv:
  contour,f,time,yy,/nodata,position=axespos,$
    xstyle=5,ystyle=5,xtick_get=xtickv
;
; 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
;
; if keyword_set(ps) then begin
;   icolor = !p.color
;   background = 255
; endif else if keyword_set(jpg) then begin
;   icolor = !p.color
;   background = 255
; endif else begin     ; X11
;   icolor = !p.color
;   background = 255
; endelse
  icolor = !p.color
  background = 255
;
; Second contour call with data and axes:
; 12/5/07: Found that the max_value does not work with the /fill 
;          keyword (does not correctly exclude missing values,
;          e.g., T,U,V at top level). Using /cell_fill instead of
;          /fill does work, but when using /cell_fill, faint gridded 
;          lines appear on ps plots (they do not show in X display 
;          device). 

; contour,f,time,yy,position=axespos,/noerase,/cell_fill,ytitle=ytitle,$
;   title=title,xminor=xminor,charsize=chsize,xstyle=1,ystyle=1,$
;   xtickformat='LABEL_DATE',c_colors=c_colors,levels=clevels,$
;   background=background,color=icolor,min_value=fmin,max_value=fmax,$
;   /follow


; Axes, and plot color fill:
  contour,f,time,yy,position=axespos,/noerase,/cell_fill,/follow,$
    xstyle=1,ystyle=1,levels=clevels,color=icolor,xtickformat='LABEL_DATE',$
    ytitle=ytitle,title=title,xminor=xminor,charsize=chsize,c_colors=c_colors,$
    min_value=fmin,max_value=fmax

; Add contour lines and labels:
  if (contour_label) then $
    contour,f,time,yy,position=axespos,/noerase,/follow,$
      levels=clevels,background=background,color=icolor,min_value=fmin,$
      max_value=fmax,xstyle=5,ystyle=5,xtickformat='LABEL_DATE'
;
; 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])
      slt = calcslt(uthr,glon)
      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
      if keyword_set(ps) then $
        xyouts,xnorm[0,0],yll-.050,sltlab,/norm,align=0.5,charsize=chsize-0.2 $
      else $ 
        xyouts,xnorm[0,0],yll-.055,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:
;
    xyouts,xmid,yll-.28,filelabel,/norm,align=0.5,charsize=chsize
;
; Color bar:
;
    form = "(g10.4)"
    cblabels = strarr(nlevels+1)
    cblabels[*] = ' '
    cblabels[0] = strcompress(string(format=form,cmin),/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,cmax),$
      /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
;
; Report min,max,interval to stdout and external text output file:
; pro myprint,text,stdout=stdout,outfile=outfile
;
;   state = *stateptr
    text = '            Contour '+minmaxlabel
    myprint,text,/stdout,outfile=state.outfile,/blankline
end ; pro utcontour