load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

;;begin                    ; optional
  diri      = "./"
  filenames = systemfunc(" cd "+diri +"; ls Sod.*.tab")
;;filenames = (/"Sod.0000.tab","Sod.0001.tab","Sod.0002.tab"/)
  nfiles    = dimsizes(filenames)
  Nx1       = new(nfiles,integer)
  time      = new(nfiles,float)
  lev       = new(nfiles,integer)
  dom       = new(nfiles,integer)

;---Loop through files to get total Nx1 count and time/level/domain for each file
  do i=0,nfiles-1
    print("--------------------------------------------------")
    print("Reading file '" + filenames(i) + "'...")    

    lines  = asciiread(filenames(i),-1,"string")
    Nx1(i) = toint(str_get_field(lines(0),4," "))
    time(i)= tofloat(str_get_field(lines(2), 6," ,="))
    lev(i) = toint(str_get_field(lines(2)  , 8," ,="))
    dom(i) = toint(str_get_field(lines(2)  ,10," ,="))
   ;delete(lines)          ; use if lines were to change size
  end do

  print(" ")
  print(Nx1+"  "+time+"  "+lev+"  "+dom)
  print(" ")

  Nx1_total = sum(Nx1)  

;---Create arrays to hold all data at each time
  izone = new((/nfiles,Nx1_total/),float)
  x1    = new((/nfiles,Nx1_total/),float)
  d     = new((/nfiles,Nx1_total/),float)
  V1    = new((/nfiles,Nx1_total/),float)
  V2    = new((/nfiles,Nx1_total/),float)
  V3    = new((/nfiles,Nx1_total/),float)
  P     = new((/nfiles,Nx1_total/),float)

  ncol  = 7    ; # of columns
  opt   = 4    ; # of header lines to skip


;---Loop through files to read data into arrays
  ibeg = 0
  do i=0,nfiles-1

    data = readAsciiTable(filenames(i),ncol,"float",opt)

    iend               = ibeg + Nx1(i) - 1
    izone(i,ibeg:iend) = data(:,0)
    x1(i,ibeg:iend)    = data(:,1)
    d(i,ibeg:iend)     = data(:,2)
    V1(i,ibeg:iend)    = data(:,3)
    V2(i,ibeg:iend)    = data(:,4)
    V3(i,ibeg:iend)    = data(:,5)
    P(i,ibeg:iend)     = data(:,6)

;---Increment for next time in loop.
    ibeg = iend + 1
  end do

  print("Nx1           = " + Nx1)
  print("min/max izone = " + min(izone) + "/" + max(izone))
  print("min/max x1    = " + min(x1) + "/" + max(x1))
  print("min/max d     = " + min(d) + "/" + max(d))
  print("min/max V1    = " + min(V1) + "/" + max(V1))
  print("min/max V2    = " + min(V2) + "/" + max(V2))
  print("min/max V3    = " + min(V3) + "/" + max(V3))
  print("min/max P     = " + min(P) + "/" + max(P))

  wks = gsn_open_wks("x11","Sod")  ; Open an X11 workstation.

  res                        = True

  res@tiMainString           = "density profile"
  res@tiXAxisString          = "position"
  res@tiYAxisString          = "density"
  res@tiXAxisFontHeightF     = 0.02        ; Change the font size.
  res@tiYAxisFontHeightF     = 0.02

 ;res@xyLineColors           = (/3,4/)     ; Set the line colors.
 ;res@xyLineThicknessF       = 2.0         ; Double the width.
 ;res@xyLineLabelFontHeightF = 0.02        ; Font size and color
 ;res@xyLineLabelFontColor   = 2           ; for line labels

  do nf=0,nfiles-1
     res@gsnCenterString = "time="+time(nf)+" level="+lev(nf)+"  domain="+dom(nf)
     plot = gsn_csm_xy(wks,x1(nf,:),d(nf,:),res) 
  end do

  res@gsnCenterString = "All Times" 
  plot = gsn_csm_xy(wks,x1,d,res) 

;;end    ; required only if "begin" is present