;-----------------------------------------------------------------------
;
; B. Foster 9/03:
;
; Idl procedures to read and make line plots from xyut netcdf file saved
;   by f90 post-processor tgcmproc (set namelist sendcdf_xyut). 
;   Fields were saved by tgcmproc at user selected locations and levels 
;   for each model time. This code makes plots with time on x-axis and 
;   the field on the y-axis at each location and level.
;
; Main pro: pro readnc_xyut,ncfile,ps=ps
;   Where ncfile is the netcdf file to read, and ps is optional keyword.
; By default (ps keyword not set), plots are displayed on X screen. 
;   In this mode, left-click in window to display next plot. 
; To get a ps file (xyut.ps) instead of displaying on the screen, call 
;   with keyword ps set (e.g.: "readnc_xyut,'xyut.nc',/ps")
;
;-----------------------------------------------------------------------
pro varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,var,$
  shortname,longname
common pltcom, iframe,mkps
;
; Plot var at all locations and levels (assume location type is lat,lon):
; (Field variables are dimensioned (nsloc,nslev,ntime))
;
  xtitle = 'CALENDAR DAY'
  pos = [0.17,0.17, 0.87,0.87]
  for l=0,nsloc-1 do begin
    for k=0,nslev-1 do begin
      if ftype ne 'HT-INDEP' then begin
        if slev_types[k] ne 'hmf2' then begin
          title=string(format=$
            "(a,':  Lat, Lon=',f8.2,',',f8.2,' ',a,' ',a,'=',f8.2)",$
            shortname,sloc(*,l),sloc_names[l],slev_types[k],slev[k]) 
        endif else begin
          title=string(format=$
            "(a,':  Lat, Lon=',f8.2,',',f8.2,' ',a,' at hmf2 heights')",$
            shortname,sloc(*,l),sloc_names[l])
        endelse
      endif else begin ; height-independent field
        title=string(format=$
          "(a,':  Lat, Lon=',f8.2,',',f8.2,' ',a)",$
          shortname,sloc(*,l),sloc_names[l])
      endelse
      iframe = iframe+1
      if mkps eq 0 then cursor,x,y,/down
      print,format="('Frame ',i3,': ',a)",iframe,title
      plot,day,var(l,k,*),title=title,xtitle=xtitle,ytitle=longname,$
        position=pos,/norm
      if ftype eq 'HT-INDEP' then goto,skiplev
    endfor ; slev loop
    skiplev:
  endfor ; sloc loop
end
;-----------------------------------------------------------------------
function varget,ncid,sname,longname,ftype
;
; Get variable sname from file attached to ncid. Also return longname
; and ftype.
;
  ncdf_varget,ncid,ncdf_varid(ncid,sname),var
  ncdf_attget,ncid,ncdf_varid(ncid,sname),'long_name',name 
  ncdf_attget,ncid,ncdf_varid(ncid,sname),'units',units 
  longname = string(name) + ' (' + string(units) + ')'
  ncdf_attget,ncid,ncdf_varid(ncid,sname),'field_type',ftype 
  ftype = string(ftype)
  return,var
end
;-----------------------------------------------------------------------
pro readnc_xyut,ncfile,ps=ps
common pltcom, iframe,mkps
;
; Open netcdf file and print global attributes:
;
  ncid = ncdf_open(ncfile,/NOWRITE)
  print,' ' & print,'ncfile=',ncfile
  inq = ncdf_inquire(ncid)
  print,'Global Attributes:'
  for i=0,inq.ngatts-1 do begin
    attname = ncdf_attname(ncid,/GLOBAL,i)
    attinq = ncdf_attinq(ncid,/GLOBAL,attname)
    ncdf_attget,ncid,/GLOBAL,attname,attval
    if attinq.datatype eq 'CHAR' then $
      print,'  ',attname,' = ',string(attval) $
    else print,'  ',attname,' = ',attval
  endfor
;
; Get dimensions:
;
; ntime = number of times (unlimited dimension)
;
  id_ntime = ncdf_dimid(ncid,'time')
  ntime_name = "" & ncdf_diminq,ncid,id_ntime,ntime_name,ntime
  print,'ntime=',ntime
;
; nslev = number of selected levels:
;
  id_nslev = ncdf_dimid(ncid,'nslev')
  nslev_name = "" & ncdf_diminq,ncid,id_nslev,nslev_name,nslev
  print,'nslev=',nslev
;
; nsloc = number of selected locations:
;
  id_nsloc = ncdf_dimid(ncid,'nsloc')
  nsloc_name = "" & ncdf_diminq,ncid,id_nsloc,nsloc_name,nsloc
  print,'nsloc=',nsloc
;
; slev = selected levels (floats):
;
  slev = fltarr(nslev)
  ncdf_varget,ncid,ncdf_varid(ncid,'slev'),slev
;
; slev_type = types of selected levels (strings):
;
  levtypes = strarr(nslev)
  ncdf_varget,ncid,ncdf_varid(ncid,'slev_type'),levtypes
  slev_types = strarr(nslev)
  for i=0,nslev-1 do slev_types[i] = string(levtypes[*,i])
;
; sloc = selected locations (float(2,nsloc))
;
  sloc = fltarr(2,nsloc)
  ncdf_varget,ncid,ncdf_varid(ncid,'sloc'),sloc
;
; sloc_type = types of selected locations (strings):
;
  loctypes = strarr(nsloc)
  ncdf_varget,ncid,ncdf_varid(ncid,'sloc_type'),loctypes
  sloc_types = strarr(nsloc)
  for i=0,nsloc-1 do sloc_types[i] = string(loctypes[*,i])
;
; sloc_names = names of selected locations:
;
  locnames = strarr(nsloc)
  ncdf_varget,ncid,ncdf_varid(ncid,'sloc_name'),locnames
  sloc_names = strarr(nsloc)
  for i=0,nsloc-1 do sloc_names[i] = string(locnames[*,i])
;
; mtimes = model times (3,ntime) (1st dimension is day,hour,minute):
; Construct decimal days for x-axis from mtimes.
;
  mtimes = intarr(3,ntime) 
  ncdf_varget,ncid,ncdf_varid(ncid,'mtime'),mtimes
  day = fltarr(ntime)
  for i=0,ntime-1 do $
    day[i] = float(mtimes(0,i))+float(mtimes(1,i))/24.+float(mtimes(2,i))/(24.*60.)
  iframe = 0
;
; Init ps device if necessary:
;
  mkps = 0
  if keyword_set(ps) then begin
    psfile = 'xyut.ps'
    set_plot,'ps'
    device,filename=psfile
    mkps = 1
  endif else set_plot,'X'
;
; Get and plot desired fields:
;
  tn = varget(ncid,'TN',tn_name,ftype) 
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,tn,'TN',tn_name
;
  un = varget(ncid,'UN',un_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,un,'UN',un_name
;
  vn = varget(ncid,'VN',vn_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,vn,'VN',vn_name
;
  w = varget(ncid,'W',w_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,w,'W',w_name
;
  o2 = varget(ncid,'O2',o2_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,o2,'O2',o2_name
;
  o1 = varget(ncid,'O1',o1_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,o1,'O1',o1_name
;
  n2 = varget(ncid,'N2',n2_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,n2,'N2',n2_name
;
  ui = varget(ncid,'UI',ui_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,ui,'UI',ui_name
;
  vi = varget(ncid,'VI',vi_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,vi,'VI',vi_name
;
  wi = varget(ncid,'WI',wi_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,wi,'WI',wi_name
;
  nelec = varget(ncid,'NE',ne_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,nelec,'NE',ne_name
;
  hmf2 = varget(ncid,'HMF2',hmf2_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,hmf2,'HMF2',hmf2_name
;
  fof2 = varget(ncid,'FOF2',fof2_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,fof2,'FOF2',fof2_name
;
  z = varget(ncid,'Z',z_name,ftype)
  varplt,day,nsloc,sloc,sloc_names,nslev,slev,slev_types,ftype,z,'Z',z_name
;
  if keyword_set(ps) then begin
    device,/close_file
    print,'Created postscript file ',psfile
  endif
end