;-----------------------------------------------------------------------
pro plot_3d,ncid,file,f,fname,dnames,dsizes,dimids,psfile
  print,format="(/,'Plot 3d var ',a,' (',2(a,'=',i3,','),a,'=',i3,')')", $
    fname,dnames[0],dsizes[0],dnames[1],dsizes[1],dnames[2],dsizes[2]

  idv = ncdf_varid(ncid,fname)
  ncdf_attget,ncid,idv,'units',units
  units = string(units)
  if (units eq 'cm') then begin
    f = f*1.e-5
    units = 'km'
    print,'Converted ',fname,' from cm to km'
  endif

  idv = ncdf_varid(ncid,dnames[0])
  ncdf_varget,ncid,idv,x
  ncdf_attget,ncid,idv,'long_name',xname
  xname = string(xname)

  idv = ncdf_varid(ncid,dnames[1])
  ncdf_varget,ncid,idv,y
  ncdf_attget,ncid,idv,'long_name',yname
  yname = string(yname)

  idv = ncdf_varid(ncid,dnames[2])
  ncdf_varget,ncid,idv,time
  ntime = n_elements(time)

 ;x0=.1 & y0=.15 & x1=.93 & y1=.93
  x0=.1 & y0=.25 & x1=.93 & y1=.93
  position = [x0,y0,x1,y1]
  ans = ''
  charsize = 1.5
  c_charsize = 1.4
  if (psfile ne '') then begin
    charsize = 1.0
    c_charsize = 1.0
  endif
  for n=0,ntime-1 do begin
    prompt = string(format="('Plot ',a,' at itime ',i4,' (of ',i4,')? (y/n/q): ')",$
      fname,n+1,ntime)
    read,ans,prompt=prompt
    case ans of
      'y': 
      'n': goto,endloop
      'q': goto,quitloop
      else: 
    endcase
    fplot = f[*,*,n]
    title = ' Field: ' + fname + ' (' + units + ') Time=' + string(time[n])
    fmin = min(fplot) & fmax = max(fplot)
    minmax_label = 'Min,Max = ' + string(fmin) + string(fmax)
    file_label = 'File = ' + file
    print,title,' min,max=',fmin,fmax
    contour,fplot,x,y,title=title,xtitle=xname,ytitle=yname,position=position,/follow,$
      charsize=charsize,c_charsize=c_charsize,nlevels=12,xstyle=1,ystyle=1
    xyouts,0.5,y0-.15,minmax_label,/normal,alignment=0.5,charsize=charsize
    xyouts,0.5,y0-.20,file_label,/normal,alignment=0.5,charsize=charsize
;   cursor,xc,yc,wait=1
    endloop:
  endfor ; n=0,ntime-1
  quitloop:
end
;-----------------------------------------------------------------------
pro plot_4d,ncid,file,f,fname,dnames,dsizes,dimids,klevs,psfile
  print,format="(/,'Plot 4d Var ',a,' (',3(a,'=',i3,','),a,'=',i3,')')", $
    fname,dnames[0],dsizes[0],dnames[1],dsizes[1],dnames[2],dsizes[2],$
    dnames[3],dsizes[3]

  idv = ncdf_varid(ncid,fname)
  ncdf_attget,ncid,idv,'units',units
  units = string(units)
  if (units eq 'cm') then begin
    f = f*1.e-5
    units = 'km'
    print,'Converted ',fname,' from cm to km'
  endif
  
  idv = ncdf_varid(ncid,dnames[0])
  ncdf_varget,ncid,idv,x
  ncdf_attget,ncid,idv,'long_name',xname
  xname = string(xname)

  idv = ncdf_varid(ncid,dnames[1])
  ncdf_varget,ncid,idv,y
  ncdf_attget,ncid,idv,'long_name',yname
  yname = string(yname)

  idv = ncdf_varid(ncid,dnames[2])
  ncdf_varget,ncid,idv,z
  ncdf_attget,ncid,idv,'long_name',zname
  zname = string(zname)
  nz = n_elements(z)

  idv = ncdf_varid(ncid,dnames[3])
  ncdf_varget,ncid,idv,time
  ntime = n_elements(time)

  nklevs = n_elements(klevs)

; x0=.1 & y0=.15 & x1=.93 & y1=.93
  x0=.1 & y0=.25 & x1=.93 & y1=.93
  position = [x0,y0,x1,y1]
  ans = ''
  charsize = 1.5
  c_charsize = 1.4
  if (psfile ne '') then begin
    charsize = 1.0
    c_charsize = 1.0
  endif
;
; Time loop:
  for n=0,ntime-1 do begin
    prompt = string(format="('Plot ',a,' at itime ',i4,' (of ',i4,')? (y/n/q): ')",$
      fname,n+1,ntime)
    read,ans,prompt=prompt
    case ans of
      'y': 
      'n': goto,endtloop
      'q': goto,quitloop
      else: 
    endcase
;
; Requested levels loop:
    for k=0,nklevs-1 do begin
      if (klevs[k] gt nz-1) then begin
        print,'Note: skipping k-level ',klevs[k],' because it is greater than nlev=',nz
        continue
      endif
      prompt = string(format="('Plot ',a,' at k-level ',i4,' (of ',i4,')? (y/n/q): ')",$
        fname,klevs[k],nz)
      read,ans,prompt=prompt
      case ans of
        'y': 
        'n': goto,endkloop
        'q': goto,endtloop
        else: 
      endcase
      fplot = f[*,*,klevs[k],n]
      title = ' Field: ' + fname + ' (' + units + ') Level=' + string(z[klevs[k]]) + ' Time=' + string(time[n])
      fmin = min(fplot) & fmax = max(fplot)
      minmax_label = 'Min,Max = ' + string(fmin) + string(fmax)
      file_label = 'File = ' + file
      print,title,' min,max=',fmin,fmax
      contour,fplot,x,y,title=title,xtitle=xname,ytitle=yname,position=position,/follow,charsize=charsize,$
        c_charsize=c_charsize,nlevels=12,xstyle=1,ystyle=1
      xyouts,0.5,y0-.15,minmax_label,/normal,alignment=0.5,charsize=charsize
      xyouts,0.5,y0-.20,file_label,/normal,alignment=0.5,charsize=charsize
;     cursor,xc,yc,wait=1
      endkloop:
    endfor ; k=0,nklevs-1
    endtloop:
  endfor ; n=0,ntime-1
  quitloop:
end
;-----------------------------------------------------------------------
pro ncplot,file=file,psfile=psfile
;
; Path to input file:
;
  if (keyword_set(file)) then begin
    file_path = file
  endif else begin
    file_path = 'waccm_edyn.nc'
;   file_path = 'timegcm_edyn.nc'
  endelse

  if (keyword_set(psfile)) then begin
    set_plot,'PS'
    device,filename=psfile
    print,'Will send ps plots to ',psfile
  endif else begin
    psfile = ''
    set_plot,'X'
    device,retain=2
  endelse
;
; Plot these 3d or 4d fields:
; (if 'all', then plot all 3d and 4d fields)
;
; fnames = ['Z','PED_MAG','trash','SINI_MAG']
; fnames = ['Z_DYN','ZPOT_MAG']
  fnames = ['all']
;
; Plot 4d fields at these k-levels:
;
  klevs = [5,15,25,35,45,55,65,75] ; k-levels to plot 

  ncid = ncdf_open(file_path) ; default is nowrite
  print,'Opened file ',file_path

  nf = n_elements(fnames)
  file = ncdf_inquire(ncid)
  name = '' & ntime=0
  for i=0,nf-1 do begin
    found = 0
    for n=0,file.nvars-1 do begin
      var = ncdf_varinq(ncid,n)
      if (var.ndims lt 2) then begin
;       print,'>>> Skipping var ',var.name,' because it has only ',var.ndims,' dimensions'
        continue
      endif 
      if (var.name eq fnames[i] or fnames[i] eq 'all') then begin ; found var on file
        if (var.name eq 'mtime') then continue
        found = 1
;       print,'Plotting var ',var.name
        ncdf_varget,ncid,n,fvar
        dimnames = strarr(var.ndims)
        dimsizes = intarr(var.ndims)
        for idims=0,var.ndims-1 do begin
          ncdf_diminq,ncid,var.dim[idims],name,size
          dimnames[idims] = name
          dimsizes[idims] = size 
        endfor 
        if (var.ndims eq 3) then begin
          plot_3d,ncid,file_path,fvar,var.name,dimnames,dimsizes,var.dim,psfile
        endif else if(var.ndims eq 4) then begin
          plot_4d,ncid,file_path,fvar,var.name,dimnames,dimsizes,var.dim,klevs,psfile
        endif
      endif ; is a requested var
    endfor ; n=0,file.nvars-1 (vars on the file)
    if (not found) then print,'>>> Could not find variable ',fnames[i]
  endfor ; i=0,nf-1 (requested vars)

  ncdf_close,ncid
  if (keyword_set(psfile)) then begin
    device,/close_file
    print,'Closed ps file ',psfile
  endif
end
;-----------------------------------------------------------------------