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/shea_util.ncl"

begin
;************************************************
; Variable and file handling
; The is available from: 
;     http://www.ncl.ucar.edu/Applications/Data/
;************************************************
   fn  = "WCA_areaaveraged.nc"                       ; define filename
   in  = addfile(fn,"r")                    ; open netcdf file
   precip  = in->PREC(:,0,0)                              ; get data

d   = 1    ; detrending opt: 0=>remove mean 1=>remove mean + detrend
sm  = 21   ; smooth: should be at least 3 and odd
pct = 0.1 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common. 

spec = specx_anal(precip,d,sm,pct)

;************************************************
; calculate confidence interval [here 5 and 95%]
; return 4 curves to be plotted
;************************************************

splt = specx_ci(spec, 0.05, 0.95)
printVarSummary(splt)                          ; splt(4,702)


;************************************************
; plotting
;************************************************
   wks  = gsn_open_wks("png","WCA")              ; Opens a ps file

   r = True                                      ; plot mods desired
   r@tiMainString = "WCA"                ; title
   r@tiXAxisString = "Frequency"  ; xaxis
   r@tiYAxisString = "Variance/frq_interval"     ; yaxis
   r@xyLineColors        = (/"foreground","green","blue","red"/)
  
   f = spec@frq                                 ; convenience
   ;p = 1/f
   ;p!0 = "f"
   ;p&f = f
   ;p@long_name = "period"
   ;p@units     = "days"
   ;r@trYMinF             = 0.00                 ; manually set lower limit
   ;r@trYMaxF             = 20.0   
   ;ip   = ind(p.le.90) 
   
   plot = gsn_csm_xy(wks,f,splt,r)     ; create plot
   ;r@trXReverse = True
   ;plot = gsn_csm_xy(wks,p(ip), splt(:,ip),r)  
end