;;   Compilation scattplot of weekly totals for every location
begin

;--- read all scatter data files
   file_path       = "/wind01/zrieck/CloudSeeding/precip/Compilations/"
   f_list = systemfunc("ls "+file_path + "Data*")
   f_in = addfiles(f_list,"r")

;  Import stations
   AC_wrf = f_in[0]->AC_wrf_prec
   AC_sn = f_in[0]->AC_Snotel_prec

   BR_wrf = f_in[1]->BR_wrf_prec
   BR_sn = f_in[1]->BR_Snotel_prec

   HR_wrf = f_in[2]->HR_wrf_prec
   HR_sn = f_in[2]->HR_Snotel_prec

   PC_wrf = f_in[3]->PC_wrf_prec
   PC_sn = f_in[3]->PC_Snotel_prec

   TL_wrf = f_in[4]->wrf_prec
   TL_sn = f_in[4]->Snotel_prec

;  Create combo array for plot- 1 for wrf and 1 for Snotel values
   wrf_Totals = new((/5,dimsizes(AC_wrf)/),float)
   sn_Totals = new((/5,dimsizes(AC_sn)/),float)
   wrf_Totals(0,:) = AC_wrf
   sn_Totals(0,:) = AC_sn
   wrf_Totals(1,:) = BR_wrf
   sn_Totals(1,:) = BR_sn
   wrf_Totals(2,:) = HR_wrf
   sn_Totals(2,:) = HR_sn
   wrf_Totals(3,:) = PC_wrf
   sn_Totals(3,:) = PC_sn
   wrf_Totals(4,:) = TL_wrf
   sn_Totals(4,:) = TL_sn

;  String into 1-D and calculate regression coefficient
   wrf_totals = ndtooned(wrf_Totals)
   sn_totals = ndtooned(sn_Totals)
   rc = regline(wrf_totals,sn_totals)
   print(rc)

   reg_data = new((/2,dimsizes(wrf_totals)/),float)
   reg_data(0,:) = wrf_totals
   reg_data(1,:) = rc * sn_totals

;  Create plot
   wks = gsn_open_wks("pdf","Combined_Scatter2")
   res = True
   res@tiMainString = "Combined Scaterplot"
   res@gsnMaximize         = True
   res@xyMarkLineModes     = (/"Markers","Lines"/)
   res@xyMarkers           = 16
   res@xyMarkerColor       = "red"
   res@xyMarkerSizeF       = 0.005
   res@xyDashPatterns      = 1             
   res@xyLineThicknesses   = (/1,2/)      
   res@tiXAxisString = "Snotel (mm)"
   res@tiYAxisString = "WRF (mm)"
   plot1 = gsn_csm_xy(wks,sn_totals,reg_data,res)

exit   
end