;
; skewt_dns.ncl: ncl routines for plotting average dropsonde data on a skew-T
;
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "~/Library/NCL/skewt_func.ncl"
;**************************************************
 begin

   source = "rvsj"
   fi = addfile("sondes_"+source+"_3.0.nc","r")
;
   rico_day = fi->rico_day
;
; define coordinates of uniform grid and a template variable for gridded data
;
   rday = fspan( -24.5,24.5,150)
   rday@long_name = "Days after 2005-01-01"
   rday@units = "days"

   plvl = new( (/dimsizes(fi->DeltaP)/), "float")
   plvl = (/avg(fi->SurfacePressure) -fi->DeltaP/)
   plvl@long_name = "Pressure"
   plvl@units = "Pa"

   template = new ( (/dimsizes(plvl),dimsizes(rday)/), "float")
   template@_FillValue = -999.
   template!0 = "plvl"
   template&plvl = plvl
   template!1 = "rday"
   template&rday = rday
;
; define coordinates of input data 
;
   xi   = rico_day
   yi   = avg(fi->SurfacePressure) - fi->DeltaP
   yi@long_name = "Pressure"
   xi@long_name = rday@long_name
   yi@units = "Pa"

   fldi = (/fi->tair/)
   fldi@_FillValue = -999.
   fldi!0 = "yi"
   fldi&yi = yi
   fldi!1 = "xi"
   fldi&xi = xi
;
; construct a pressure field for plvl data
;
   prsi = (/fi->tair/)
   do k=0,dimsizes(xi)-1
     prsi(:,k) = fi->SurfacePressure(k) - fi->DeltaP
   end do
;
; create plots
;
; order: Surface is 'bottom'  eg: 1000,950,935,897,...  
                          
  npts = dimsizes(prsi(:,0)) - 8
  p    = dim_avg(prsi(0:npts,:))/100.
  tc   = dim_avg(fi->tair(0:npts,:)) - 273.15
  tdc  = dim_avg(dewtemp_trh(fi->tair(0:npts,:),fi->rhum(0:npts,:))) - 273.15  ; dew pt temp  [C]   
  z    = dim_avg(fi->galt(0:npts,:))    ; geopotential [gpm] 
  wspd = dim_avg(fi->wspd(0:npts,:))    ; wind speed   [knots or m/s]    
  wdir = dim_avg(fi->wdir(0:npts,:))    ; meteorological wind dir   
;
;*************************
; create plot
;*************************
;
  wks  = gsn_open_wks ("pdf", "skewt_"+source)

; --- Create background skew-T and plot sounding----------------

  skewtOpts                 = True
  skewtOpts@DrawColAreaFill = True    ; default is False
  skewtOpts@tiMainString    = "ATS Rawindsonde - default dataOpts" 


  dataOpts           = True     ; options describing data and ploting
  dataOpts@Wthin     = 10        ; plot every n-th wind barb

  skewtOpts@tiMainString    = "RICO "+source+" Sondes" 
  skewtOpts@DrawFahrenheit  = False   ; default is True
  skewtOpts@DrawColAreaFill = True

  skewt_bkgd = skewT_BackGround (wks, skewtOpts)
  skewt_data = skewT_PlotData   (wks, skewt_bkgd, p,tc,tdc,z \
                                    , wspd,wdir, dataOpts)
;  draw (skewt_bkgd)
  draw (skewt_data)
  frame(wks)

 end