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

begin

; Open and read our station data. We know there are 26 columns (lat, lon, and 24 hours) and 73 stations
  nrow =  73
  ncol = 26
  amtdat = asciiread("data/StationDataColorado_amount.txt",(/nrow,ncol/),"float")
  occdat = asciiread("data/StationDataColorado_occurrence.txt",(/nrow,ncol/),"float")
  intdat = asciiread("data/StationDataColorado_intensity.txt",(/nrow,ncol/),"float") 

; Set lats and lons for each station
  lon = amtdat(:,0)
  lat = amtdat(:,1)

; Create new variables for the station data
  pamtob = new((/24,73/),float)
  poccob = new((/24,73/),float)
  pintob = new((/24,73/),float)

; Loop through station data and write to our new variables we created.
  j = 0
  do i=2,25
    pamtob(j,:) = amtdat(:,i)
    poccob(j,:) = occdat(:,i)
    pintob(j,:) = intdat(:,i)
    j = j+1
  end do

; Average over each station in Colorado
  pintob = where(pintob.eq.-999.,pintob*0,pintob)
  amtob = dim_avg_n(pamtob,1)
  occob = dim_avg_n(poccob,1)
  intob = dim_avg_n(pintob,1)

; Create new variables for the data we are plotting. We will plot the station data and data from three
;  ensemble members
  amt = new((/4,24/),float)
  occ = new((/4,24/),float)
  int = new((/4,24/),float)

; Writing the station data to the new arrays we created.
  amt(0,:) = amtob
  occ(0,:) = occob
  int(0,:) = intob

; Here are the ensemble members we want to look at
  mem = (/"rnty","rktm","rtty"/)
  numMEM = dimsizes(mem)

; Start to loop through each ensemble member to write the data to our new arrays.
  do m=0,numMEM-1
; Read in the files we made in the previous step.
    a = addfile("data/prec_diurnal_jja_"+ mem(m) +".nc","r")

; Create new arrays for the model data
    pmodp = new((/24,73/),float)
    pmodpocc = new((/24,73/),float)
    pmodpint = new((/24,73/),float)

; Reading in the model data
    modp = a->PREC
    modpocc = a->PRECOCC
    modpint = a->PRECINT

    llres = True
    llres@returnInt = True

; Here we are looping through each station to find what grid cell each station is in. We use the
;  function wrf_user_ll_to_ij which does this for us. We then take the data from that grid cell.
    do l=0,nrow-1
      locij = wrf_user_ll_to_ij(a,lon(l),lat(l),llres)
      locij = locij-1
      locX = locij(0)
      locY = locij(1)
      pmodp(:,l) = modp(:,locY,locX)
      pmodpocc(:,l) = modpocc(:,locY,locX)
      pmodpint(:,l) = modpint(:,locY,locX)
    end do

; Here we are creating to arrays to write the data to the correct time zone. Our station data is 
;  in mountain time, but our model data is UTC. Therefore, we need to shift the model data to 
;  plot our diurnal cycle correctly.
    pmodpn = new((/24,73/),float)
    pmodpoccn = new((/24,73/),float)
    pmodpintn = new((/24,73/),float)
    pmodpn(0:16,:) = pmodp(7:23,:)
    pmodpn(17:23,:) = pmodp(0:6,:)
    pmodpoccn(0:16,:) = pmodpocc(7:23,:)
    pmodpoccn(17:23,:) = pmodpocc(0:6,:)
    pmodpintn(0:16,:) = pmodpint(7:23,:)
    pmodpintn(17:23,:) = pmodpint(0:6,:)

; Here we average over each station for the model data and calculate the intensity.
    amtmod = dim_avg_n(pmodpn,1)
    occmod = dim_avg_n(pmodpoccn,1)
    intmod = amtmod/occmod

; Write the model data to our array that includes the station data
    amt(m+1,:) = amtmod
    occ(m+1,:) = occmod
    int(m+1,:) = intmod

    delete([/amtmod,occmod,intmod,pmodpn,pmodpoccn,pmodpintn,pmodp,pmodpocc,pmodpint/])
    delete([/modp,modpocc,modpint,a/])
  end do

  tau = new(24,integer)
  do ii=0,23
    tau(ii) = ii
  end do

; Open our workstation to make the plot
  wks = gsn_open_wks("X11","precip_diurnal_xy_panel")
  plot = new(3,graphic)
; Choosing our options for the time series plot
  res = True
  res@gsnDraw = False
  res@gsnFrame = False
  res@tmXTOn            = False         ; turn off the top tick marks
  res@tmYROn            = False
  res@tiXAxisString = "Hour"
  res@trXMaxF = 23.
  res@xyLineColors = (/"black","red","blue","green"/)
  res@xyDashPatterns = (/0,0,0,0/)
  res@xyLineThicknesses = (/2.0,1.0,1.0,1.0/)

; Here we specify specifics for the legend
  res@pmLegendDisplayMode    = "Always"
  res@xyExplicitLegendLabels = (/"OBS","RNTY","RKTM","RTTY"/)
  res@pmLegendSide           = "Top"               ; Change location of
  res@pmLegendParallelPosF   = .895                ; move units right
  res@pmLegendOrthogonalPosF = -0.14               ; 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

; Here we are making Y axis labels for each panel.
  res1 = res
  res1@tiYAxisString = "Amount (mm)"
  res2 = res
  res2@tiYAxisString = "Occurences"
  res3 = res
  res3@tiYAxisString = "Intensity (mm/hour)"

; Here we make our xy plots for amount, occurence, and intensity.
  plot(0) = gsn_xy(wks,tau,amt,res1)
  plot(1) = gsn_xy(wks,tau,occ,res2)
  plot(2) = gsn_xy(wks,tau,int,res3)

; Here we plot the panels and the main title.
  resP = True
  resP@gsnFrame = False
  resP@txString = "Diurnal Cycle of Precipitation over Colorado - JJA"
  gsn_panel(wks,plot,(/1,3/),resP)

  frame(wks)

end