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

begin
fili = "surf_area.nc"

f = addfile(fili,"r")

;Set variable of interest

data = f->surf_area
time = f->T

;fix the mismatched _FillValue to properly recogonize the missing values after they are read in
data@_FillValue = totype(data@missing_value,typeof(data))
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(f->T,0)     ; returns N x 6 array
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     = time(ind_subset)

wks = gsn_open_wks("png","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