;-----------------------------------------------------------------------
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
;-----------------------------------------------------------------------
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 mylatslice,f,x,y,cmin,cmax,cint,position,title,chsize
  fmin = min(f)
  fmax = max(f)
  xll=position[0] & yll=position[1] & xur=position[2] & yur=position[3]
  xmid = (xur+xll)*0.5
  clevels = mkclevels(fmin,fmax,cmin,cmax,cint)
  nlevels = n_elements(clevels)
  loadct,33,ncolors=nlevels,bottom=3,/silent
  icolor = !p.color
  background = 255
  c_colors = indgen(nlevels)+3
  xtitle = "Longitude (deg)"
  ytitle = "Height (km)"

  contour,f,x,y,/follow,/cell_fill,xstyle=1,ystyle=1,charsize=chsize,$
    title=title,xtitle=xtitle,ytitle=ytitle,levels=clevels,position=position,color=icolor,$
    background=background,c_colors=c_colors

  contour,f,x,y,/follow,xstyle=5,ystyle=5,charsize=chsize,/noerase,$
    title=title,xtitle=xtitle,ytitle=ytitle,levels=clevels,position=position,$
    background=background,c_colors=c_colors

  minmax = string(format="('Min,Max = ',e12.4,',',e12.4,' Interval=',e12.4)",$
    cmin,cmax,cint)
  xyouts,xmid,yll-.18,minmax,/norm,align=0.5

;   position = [xll,yll,xur,yur]
  barw = xur-xll-.05*(xur-xll) & barh = .03
  xbar = xmid-.5*barw
  ybar = yll-barh-.08
  baroff = .035

  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

end ; mylatslice
;-----------------------------------------------------------------------
pro nc_latslice,file=ncfile,psfile=psfile
; if (not keyword_set(file)) then ncfile  = "/aim/d/tgcm/data/gwf_oct01_1822_to_2408_UT.nc"
; if (not keyword_set(file)) then ncfile  = "/oreo/d/foster/mkgwf/gwf_oct01_1822_to_2408_UT.nc"
  ncid = ncdf_open(ncfile)
  print,'Opened nc file ',ncfile

;dimensions:
;	time = UNLIMITED ; // (54 currently)
;	lon = 111 ;
;	lat = 55 ;
;	lev = 51 ;
;
variables:
;	int day(time) ;
;	float ut(time) ;
;	float dday(time) ;
;	float time(time) ;
;	float lon(lon) ;
;	float lat(lat) ;
;	float lev(lev) ;
;
  id = ncdf_dimid(ncid,'lon')
  ncdf_diminq,ncid,id,name,nlon
  id = ncdf_varid(ncid,'lon')
  ncdf_varget,ncid,id,lon
; print,format="(/,'nlon=',i4,' lon=',/,(8f8.2))",nlon,lon

  id = ncdf_dimid(ncid,'lat')
  ncdf_diminq,ncid,id,name,nlat
  id = ncdf_varid(ncid,'lat')
  ncdf_varget,ncid,id,lat
; print,format="(/,'nlat=',i4,' lat=',/,(8f8.2))",nlat,lat

  id = ncdf_dimid(ncid,'lev')
  ncdf_diminq,ncid,id,name,nlev
  id = ncdf_varid(ncid,'lev')
  ncdf_varget,ncid,id,lev
; print,format="(/,'nlev=',i4,' lev=',/,(8f8.2))",nlev,lev

  id = ncdf_dimid(ncid,'time')
  ncdf_diminq,ncid,id,name,ntime
  id = ncdf_varid(ncid,'time')
  ncdf_varget,ncid,id,time
; print,format="(/,'ntime=',i4,' time=',/,(8f8.2))",ntime,time

  id = ncdf_dimid(ncid,'time')
  ncdf_diminq,ncid,id,name,ntime
  id = ncdf_varid(ncid,'time')
  ncdf_varget,ncid,id,time
; print,format="(/,'ntime=',i4,' time=',/,(8f8.2))",ntime,time

  id = ncdf_varid(ncid,'day')
  ncdf_varget,ncid,id,doy
; print,format="(/,'ntime=',i4,' doy=',/,(8f8.2))",ntime,doy

  id = ncdf_varid(ncid,'ut')
  ncdf_varget,ncid,id,ut
; print,format="(/,'ntime=',i4,' ut=',/,(8f8.2))",ntime,ut

  id = ncdf_varid(ncid,'dday')
  ncdf_varget,ncid,id,dday
; print,format="(/,'ntime=',i4,' dday=',/,(8f8.2))",ntime,dday

  id = ncdf_varid(ncid,'FU')
  ncdf_varget,ncid,id,fu
; print,format="(/,'fu min,max=',2f10.4)",min(fu),max(fu)

  id = ncdf_varid(ncid,'FV')
  ncdf_varget,ncid,id,fv
; print,format="('fv min,max=',2f10.4)",min(fv),max(fv)
;
; Set device and plot position:
;
  if (keyword_set(psfile)) then begin
    set_plot,'ps'
    page_width = 8.5
    page_height = 11.0
    ysize = 6.0
;   if (not keyword_set(xsize)) then xsize = 10.
    if (not keyword_set(xsize)) then xsize = 8.
    xoffset = (page_width-xsize)*0.5
    yoffset = (page_height-ysize)*0.5
    device,file=psfile,xsize=xsize,ysize=ysize,/inches,$
      xoffset=xoffset,yoffset=yoffset,/color
    xll = 0.15 & yll = 0.20
    xur = 0.90 & yur = 0.93
    position = [xll,yll,xur,yur]
    chsize = 0.8
  endif else begin
    set_plot,'X'
    xll = 0.15 & yll = 0.20
    xur = 0.90 & yur = 0.93
    position = [xll,yll,xur,yur]
    chsize = 1.2
  endelse
;
  iframe = 0
  itime0 = 20
  itime1 = 50
  delta_time = 10
  delta_lev = 2
 
; for itime=itime0,itime1,delta_time do begin
;   print,' '
;   allmin=1.e36
;   allmax=-1.e36
;   for k=0,nlev-1,delta_lev do begin 
;     iframe = iframe+1
;     title = 'FU (zonal forcing) '+$
;             ' Day = '+string(format="(i3)",doy(itime))+$
;             ' Ut = ' +string(format="(f5.2)",ut(itime))+$
;             ' Height = '+string(format="(f8.2)",lev(k))
;     fmin = min(fu[*,*,k,itime])
;     fmax = max(fu[*,*,k,itime])
;     if (fmin lt allmin) then allmin=fmin
;     if (fmax gt allmax) then allmax=fmax
;     print,format=$
;       "('Contouring FU: iframe=',i4,' day=',i3,' ut=',f5.2,' lev=',f6.2,' fu min,max=',e12.4,',',e12.4)",$
;       iframe,doy(itime),ut(itime),lev(k),fmin,fmax

;     cmin = 0. & cmax = 0. & cint = 0.
;     mycontour,fu[*,*,k,itime],lon,lat,cmin,cmax,cint,position,title,chsize
;     if (not keyword_set(psfile)) then cursor,x,y,/wait,/up

;   endfor
;   print,'dday=',dday(itime),' min,max=',allmin,allmax
; endfor

; for itime=itime0,itime1,delta_time do begin
;   print,' '
;   allmin=1.e36
;   allmax=-1.e36
;   for k=0,nlev-1,delta_lev do begin 
;     iframe = iframe+1
;     title = 'FV (zonal forcing) '+$
;             ' Day = '+string(format="(i3)",doy(itime))+$
;             ' Ut = ' +string(format="(f5.2)",ut(itime))+$
;             ' Height = '+string(format="(f8.2)",lev(k))
;     fmin = min(fv[*,*,k,itime])
;     fmax = max(fv[*,*,k,itime])
;     if (fmin lt allmin) then allmin=fmin
;     if (fmax gt allmax) then allmax=fmax
;     print,format=$
;       "('Contouring FV: iframe=',i4,' day=',i3,' ut=',f5.2,' lev=',f6.2,' fu min,max=',e12.4,',',e12.4)",$
;       iframe,doy(itime),ut(itime),lev(k),fmin,fmax

;     cmin = 0. & cmax = 0. & cint = 0.
;     mycontour,fv[*,*,k,itime],lon,lat,cmin,cmax,cint,position,title,chsize
;     if (not keyword_set(psfile)) then cursor,x,y,/wait,/up

;   endfor
;   print,'dday=',dday(itime),' min,max=',allmin,allmax
; endfor

;
; Make lat slice for test:
;
  itime = 9  ; day 274, ut 19.5
  jlat = [5,19,29,40] ; -20.0, -13.64, -9.0, -4.1

  for j=0,n_elements(jlat)-1 do begin
;   gwsfu = fltarr(nlon,nlev)
;   gwsfu[*,*] = fu[*,jlat[j],*,itime]
    gwsfu = fltarr(78,nlev) ; lon -80 to -45
    gwsfu[*,*] = fu[0:77,jlat[j],*,itime]
    cmin = 0. & cmax = 0. & cint = 0.
    title = 'FU (zonal forcing) '+$
            ' Day = '+string(format="(i3)",doy(itime))+$
            ' Ut = ' +string(format="(f5.2)",ut(itime))+$
            ' Lat = '+string(format="(f8.2)",lat(jlat[j]))
;pro mylatslice,f,x,y,cmin,cmax,cint,position,title,chsize
    print,'gwsfu: j=',j,' jlat[j]=',jlat[j],' lat=',lat[jlat[j]]
    mylatslice,gwsfu,lon[0:77],lev,cmin,cmax,cint,position,title,chsize

    gwsfv = fltarr(78,nlev) ; lon -80 to -45
    gwsfv[*,*] = fv[0:77,jlat[j],*,itime]
    cmin = 0. & cmax = 0. & cint = 0.
    title = 'FV (meridional forcing) '+$
            ' Day = '+string(format="(i3)",doy(itime))+$
            ' Ut = ' +string(format="(f5.2)",ut(itime))+$
            ' Lat = '+string(format="(f8.2)",lat(jlat[j]))
    print,'gwsfv: j=',j,' jlat[j]=',jlat[j],' lat=',lat[jlat[j]]
    mylatslice,gwsfv,lon[0:77],lev,cmin,cmax,cint,position,title,chsize
  endfor 

  if (keyword_set(psfile)) then device,/close_file
end