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


begin

;Go to directory of interest 

diri = "/scratch/conte/s/smit1770/RRB_Test_Simulation/RedRiver/OUTPUT_TWI.1hr_180402/SURFACEMODEL/2009*/"


;list names of all of the files and add ful file path
fils = systemfunc("ls "+diri+"*.d01.nc") 
  
;print(fils)   ; make sure the desired files are listed in the correct order


;read in all files  
 a = addfiles(fils,"r")
ListSetType (a, "join")
  data = a[:]->vic_lake_area_inst
  data@time=a[:]->time
  ;print(data@time)


avg_data = wgt_areaave_Wrap(data,1.,1.,0)     ; _Wrap version of function retains metadata

printVarSummary(avg_data)
printMinMax(avg_data,0)

;---Convert time to year, month, day so we can pick out desired month
date       = cd_calendar(a[:]->time,0)     ; returns N x 6 array
print(date)
ind_subset = ind(date(:,1).eq.3)     ; Grab indexes where the month is equal to March

if(ismissing(ind_subset(0))) then
  print("Didn't find value any values for March. Quitting")
  exit
end if

;---Use indexes to subset data
avg_data_subset = avg_data(ind_subset)
time_subset     = data@time(ind_subset)

wks = gsn_open_wks("png","LIS_VIC_timeseries")
 
;set resources

resplot = True

resplot@tiMainString  = "Lake Area (m^2)";title name
resplot@trXMinF       = min(time_subset)
resplot@trXMaxF       = max(time_subset)
restick = True
restick@ttFormat = "%d"
time_axis_labels(time_subset,resplot,restick)

plot = gsn_csm_xy(wks,time_subset,avg_data_subset,resplot)

end