load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"

begin
;
; Read variable of interest and fix the mismatched _FillValue to 
; properly recogonize the missing values after they are read in.
;
  fili = "VIC_SWE.nc"
  f = addfile(fili,"r")
  vic_swe = f->SWE

  vic_swe@_FillValue := totype(vic_swe@missing_value,typeof(vic_swe))
  printVarSummary(vic_swe)

;---Now read in LIS files
  diri = "./"
  lis_files = systemfunc("ls "+diri+"*.d01.nc") 

;---Parse filenames to get time and convert it to a "days since" format.
  ntime = dimsizes(lis_files)
  years = toint(str_get_cols(lis_files,11,14))
  mons  = toint(str_get_cols(lis_files,15,16))
  days  = toint(str_get_cols(lis_files,17,18))
  hours = new(ntime,integer)
  mins  = new(ntime,integer)
  secs  = new(ntime,double)
  hours = 0
  mins  = 0
  secs  = 0
;---Convert LIS time to same units as VIC time
  time = cd_inv_calendar(years,mons,days,hours,mins,secs,vic_swe&T@units,0)  

;---Read in all LIS files  
   a = addfiles(lis_files,"r")
  ListSetType (a, "join")
  lis_swe = a[:]->SWE_inst
  lis_lat = a[0]->lat
  lis_lon = a[0]->lon

;---Attach new time array to lis_swe
  lis_swe!0 = "time"
  lis_swe&time = time
  printVarSummary(lis_swe)    ; examine data array. Result is [ncl_join | 31] x [north_south | 32] x [east_west | 45]


;---Calculate average trying two different functions
;  vic_swe_avg = wgt_areaave_Wrap(vic_swe,1.,1.,0)
;  lis_swe_avg = wgt_areaave_Wrap(lis_swe,1.,1.,0)

; dim_avg_n produces the same values
  vic_swe_avg = dim_avg_n_Wrap(vic_swe,(/1,2/))
  lis_swe_avg = dim_avg_n_Wrap(lis_swe,(/1,2/))

  wks = gsn_open_wks("png","timeseries")

  res                   = True
  res@gsnMaximize       = True
  res@vpWidthF          = 0.8
  res@vpHeightF         = 0.3
  res@xyLineThicknessF  = 5.0
  res@xyMarkLineMode    = "MarkLines"
  res@xyMarker          = 16     ; filled dot

;---set resources
  res_vic = res
  res_lis = res
  restick_vic = True
  restick_lis = True
  restick_vic@ttmFormat = "%n/%d/%Y"
  restick_lis@ttmFormat = "%n/%d/%Y"
  time_axis_labels(vic_swe&T,res_vic,restick_vic)
  time_axis_labels(lis_swe&time,res_lis,restick_lis)

;
; Create two plots showing the two different average methods, which
; should be the same since we used a weights of 1.0
;
  res_vic@trXMinF           = min(vic_swe&T)
  res_vic@trXMaxF           = max(vic_swe&T)
  res_vic@xyLineColor       = "blue"
  res_vic@xyMarkerColor     = res_vic@xyLineColor

;---Plot VIC data using a lat/lon slice at 45.4375 and -96.4375
  res_vic@gsnLeftString          = "VIC: SWE"
  res_vic@gsnLeftStringFontColor = res_vic@xyLineColor
  plot_vic = gsn_csm_xy(wks,vic_swe&T,vic_swe(:,{45.4375},{-96.4375}),res_vic)

;---Compare the first three values of both arrays at what I think is the same lat/lon locations
  print("======================================================================")
  print("VIC and LIS values at lat = 45.4375 and lon = -96.4375")
  print("  VIC   /  LIS  ")
  print(sprintf("%7.4f",vic_swe(0:2,{45.4375},{-96.4375})) + " / " + \
        sprintf("%7.4f",lis_swe(:,0,25)))
  print("======================================================================")

;---Create plot of averaged VIC data
  res_vic@gsnLeftString          = "VIC: avg(SWE)"
  res_vic@gsnLeftStringFontColor = res_vic@xyLineColor
  plot_vic_avg = gsn_csm_xy(wks,vic_swe_avg&T,vic_swe_avg,res_vic)

  res_lis@trXMinF           = min(lis_swe&time)
  res_lis@trXMaxF           = max(lis_swe&time)
  res_lis@xyLineColor       = "orange"
  res_lis@xyDashPattern     = 2
  res_lis@xyMarkerColor     = res_lis@xyLineColor

;;;---CANNOT DO THIS FOR LIS DATA BECAUSE IT IS CURVILINEAR
;;  plot_lis = gsn_csm_xy(wks,lis_swe&time,lis_swe(:,{45.4375},{-96.4375}),res_lis)

;
; We had to do some debugging to figure out that index 0 & 25 of 
; lis_lat & lis_lon match with 45.4375 and -96.4375.
;
;  print(lis_lat(0,:))
;  print(lis_lon(:,25))
  res_lis@gsnRightString          = "LIS: SWE"
  res_lis@gsnRightStringFontColor = res_lis@xyLineColor
  plot_lis = gsn_csm_xy(wks,lis_swe&time,lis_swe(:,0,25),res_lis)

;---Create plot of averaged LIS data
  res_lis@gsnRightString          = "LIS: avg(SWE)"
  res_lis@gsnRightStringFontColor = res_lis@xyLineColor
  plot_lis_avg = gsn_csm_xy(wks,lis_swe_avg&time,lis_swe_avg,res_lis)
 
;---Overlay the two slice plots and draw
  overlay(plot_vic,plot_lis)
  draw(plot_vic)
  frame(wks)

;---Overlay the two averaged plots and draw
  overlay(plot_vic_avg,plot_lis_avg)
  draw(plot_vic_avg)
  frame(wks)
end