;**************************************************
; skewt_3.ncl
;
; Concepts illustrated:
;   - Drawing Skew-T plots
;   - Thinning the wind barbs in a Skew-T plot
;   - Customizing the background of a Skew-T plot
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/home/csmith/enrr/skewt_func.ncl"
;**************************************************
 begin

; --- Read Data ----------------------------------------;
  print(fili+" ")
  filein=fili
  
  titledate = fili
  head = readAsciiHead(fili, 26);
  printVarSummary(head)
  cnlvl=systemfunc("/usr/local/bin/wc "+ filein)
  nlvl=stringtointeger(cnlvl)
  nlvl=nlvl-26
    ncol = 11
    x = ispan (1,nlvl,1)
    nskip=20


if((nlvl.le.2000).and.(nlvl.ge.51))then
    nskip=30
end if

if((nlvl.le.4000).and.(nlvl.ge.2001))then
    nskip=40
end if

if((nlvl.le.8000).and.(nlvl.ge.4001))then
    nskip=60
end if

if((nlvl.ge.8001))then
    nskip=80
end if

     
  TestData = readAsciiTable(fili , (/ncol/), "float",26)
; order: Surface is 'bottom'  eg: 1000,950,935,897,...  
  TestData@_FillValue=-999.
  TestData@missing_value=-999.

    p    = TestData (:,0)    ; pressure     [mb / hPa] 
    tc   = TestData (:,1)    ; temperature  [C]
    tc=tc-273.15   
    tdc  = TestData (:,5)    ; dew pt temp  [C]
    z    = TestData (:,1)    ; geopotential [gpm] 
    wspd=p
    wspd=0.
    wdir=p
    wdir=0.
    	

;*************************
; create plot
;*************************
  wks  = gsn_open_wks ("png", outname)
  print(" "+outname)
; --- Create background skew-T and plot sounding----------------

  skewtOpts                 = True
  skewtOpts@DrawColAreaFill = True    ; default is False
  skewtOpts@tiMainString    = label
  skewtOpts@DrawWind        = False
  dataOpts                  = True    ; options describing data and ploting
  dataOpts@DrawWind         = False
  dataOpts@colWindZ ="white"
  skewt_bkgd = skewT_BackGround (wks, skewtOpts)
  i=0

  skewt_data = skewT_PlotData   (wks, skewt_bkgd, p,tc,tdc,z \
                                    , wspd,wdir, dataOpts)
  draw (skewt_bkgd)
  draw (skewt_data)
  infoTimeStamp(wks,0.015,"NOAA/ESRL PSD ~C~ Preliminary Data")
  frame(wks)


 end