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

; List the analysis we care about and read in the file to get the variable names
anal = "CFSR"
numANA = dimsizes(anal)
anaDIR = "../data"
a = addfile(anaDIR + "/" + "t2_daily_State_Data_US_cfsr.nc","r")

vNames = getfilevarnames (a) 
nNames = dimsizes (vNames)  
print(vNames)
delete(a)

; This section is taking the variable names and making them into a format we want
;   for the plot titles.  You can ignore it.
printNames = vNames
varNames = vNames
ind1 = str_index_of_substr(vNames(1) , "_", -1) + 1
ind2 = ind1-2
do ii = 1,nNames-1
  foo1 = stringtochar(vNames(ii))
  foo2 = foo1(ind1:)
  foo3 = foo1(:ind2)
  printNames(ii) = chartostring(foo2)
  varNames(ii) = chartostring(foo3)
  delete ( [/foo1,foo2,foo3/] )
end do

; The ensembles we care about and set up the directory and root of the files    
ensembles = (/"rtty_1990"/)
numENS = dimsizes(ensembles)
DIR = "../data"
rootFILE = "t2_daily_State_Data_US_"

; Number of PDFs we want to make
numPDFs = numENS + numANA

; Setting up the legend and colours of the pdf lines
legends = new(numPDFs,string)
legends = " " 
legends(0) = "RTTY_1990"
legends(1) = anal
lines = new(numPDFs,float)
lines = 2.5
colours = new(numPDFs,string)
colours = (/"blue","black"/)

; Open the workstation to make the plot
wks  = gsn_open_wks ("X11","PDF_T2_Annual_perk")
; Setting the plot resources
res = True
res@gsnDraw         = False
res@gsnFrame        = False
res@gsnMaximize     = True
res@tiYAxisString          = "PDF (%)"
res@xyLineThicknesses      = lines
res@xyLineColors           = colours
res@xyDashPatterns         = 0

res@pmLegendDisplayMode    = "Always"
res@xyExplicitLegendLabels = legends
res@pmLegendSide           = "Top"               ; Change location of
res@pmLegendParallelPosF   = .8425                 ; move units right
res@pmLegendOrthogonalPosF = -0.2325               ; move units down
res@pmLegendWidthF         = 0.125               ; Change width and
res@pmLegendHeightF        = 0.075                ; height of legend.
res@lgPerimOn              = True                ; turn off/on box around
res@lgLabelFontHeightF     = .015                ; label font height

; Set PDF resources
pdfres = True
pdfres@bin_min = -40.  ; Lowest bin is -40
pdfres@bin_max = 45.   ; Highest bin is 45
nBin = 85  ; Number of bins
; Making a variable for the plot
plots = new(nNames-1 , graphic)
txid = new(2,graphic)
amid = new(2,graphic)

; Loop through each variable for the analysis and ensembles.  Make the PDFs
do ivar =1,nNames-1
  print("Working on " + vNames(ivar) )
  do iens=0,numPDFs-1
    if (iens .ge. numENS) then
      ana = addfile(anaDIR + "/" + "t2_daily_State_Data_US_cfsr.nc","r")
      var = ana->$vNames(ivar)$
      var = var - 273.15
      delete(ana)
    else
      ensFILES = systemfunc (" ls -1 " + DIR + "/" + rootFILE + ensembles(iens) + ".nc")
      ens = addfiles(ensFILES,"r")
      var = ens[:]->$vNames(ivar)$
      delete(ensFILES)
      delete(ens)
    end if

; Make the PDF using the pdfx function
    ap  = pdfx(var, nBin, pdfres)
    if (iens .eq. 0) then
      xx = new ( (/numPDFs,nBin/), typeof(ap))
      yy = new ( (/numPDFs,nBin/), typeof(ap))
    end if
    xx(iens,:)     = ap@bin_center     ; assign appropriate "x" axis values
    yy(iens,:)     = (/ap/)
    delete(ap)
    delete(var)

  end do

; Calculate the Perkins Skill Score and print it to the screen
  t2ens_pdfx = doubletofloat(yy(0,:)/100)
  t2obs_pdfx = doubletofloat(yy(1,:)/100)
  t2_perkins = new((/nBin/),float)
  do i=0,nBin-1
    t2_perkins_temp = new((/2/),float)
    t2_perkins_temp(0) = t2ens_pdfx(i)
    t2_perkins_temp(1) = t2obs_pdfx(i)
    t2_perkins(i) = min(t2_perkins_temp)
    delete(t2_perkins_temp)
  end do
  sum_t2_perkins = sum(t2_perkins)
  print("PKSS = " +flt2string(sum_t2_perkins))

; Set the title for each panel
  res@tiMainString           = varNames(ivar) + " " +printNames(ivar)
; Make each of the PDF panels
  plots(ivar-1) = gsn_csm_xy (wks, xx, yy, res)

; Set text size and area for Perkins Skill Score to go on upper left of panels
  txres = True
  txres@txFontHeightF = 0.03
  amres = True
  amres@amParallelPosF = -0.475
  amres@amOrthogonalPosF = -0.475
  amres@amJust = "TopLeft"

; Plot Perkins Skill Sore on upper left of panels
  txid(ivar-1) = gsn_create_text(wks,"pkss = "+flt2string(sum_t2_perkins),txres)
  amid(ivar-1) = gsn_add_annotation(plots(ivar-1),txid(ivar-1),amres)
  delete( [/xx,yy/] )
end do

; Set panel plot resources
pnlres                 = True
pnlres@gsnMaximize     = True

; Draw the panel plot
gsn_panel(wks,(/plots/),(/1,2/),pnlres)

; End of program
end