;;  Brighton Code file for weekly precip total generation
begin

;--- read all CPC precipitation files
   file_path       = "/cloud01/Model/WRF/NCAR_13yr/2d_ctrl/precip/"
   f_list = systemfunc("ls "+file_path + "wrf2d_d01*")

   f_Oct = str_match_ind(f_list, "10-")
   f_Jan = str_match_ind(f_list, "01-")

   n_block = dimsizes(f_Oct)
   f_list_ind = new(n_block*2, "integer")
   ct = 0
   do i = 0, n_block-1
      f_list_ind(ct) = f_Oct(i)
      f_list_ind(ct+1) = f_Jan(i)
      ct = ct + 2
   end do

   f_in = addfiles(f_list(f_list_ind), "r")
   ListSetType(f_in, "cat")
   f_time = tostring(f_in[:]->Times)

;---use gridpoint data
;   set lat-lon coordinates to match Brighton
    a  = 40.6
    b  = -111.58

    lat2d = f_in[0]->XLAT
    lon2d = f_in[0]->XLONG

;   get x-y indices for lat-lon data, convert back to verify accuracy
    nm = getind_latlon2d(lat2d,lon2d,a,b)
    nm1 = nm(0,0)
    nm2 = nm(0,1)

;  Index values for time
   yyyy_wrf = toint(str_get_cols(f_time, 0,3))
   mm_wrf  = toint(str_get_cols(f_time, 5,6))
   dd_wrf   = toint(str_get_cols(f_time, 8,9))
   hh_wrf   = toint(str_get_cols(f_time, 11,12))

   years = ispan(2000,2013,1)
   n_years = dimsizes(years)
   mon_length = (/30,31,31,28,31/)
   n_months = dimsizes(mon_length)
   n_val = (n_months * (n_years-2))

;  Set data to skip Oct- giving Nov 1- March 31 every year
   wrf_s_mm = ind(mm_wrf.EQ.11)
   wrf_s_dd = ind(dd_wrf.EQ.01)
   wrf_s_hh = ind(hh_wrf.EQ.00)
   temp = venn2_intersection(wrf_s_dd, wrf_s_hh)
   wrf_s_ind = venn2_intersection(wrf_s_mm,temp)
   delete(temp)

   wrf_e_mm = ind(mm_wrf.EQ.03)
   wrf_e_dd = ind(dd_wrf.EQ.31)
   wrf_e_hh = ind(hh_wrf.EQ.23)
   temp = venn2_intersection(wrf_e_dd, wrf_e_hh)
   wrf_e_ind = venn2_intersection(wrf_e_mm,temp)
   delete(temp)

;  Create an array with precip data between these indices organized by row
   n_hr_season = wrf_e_ind(7) - wrf_s_ind(7) + 1
   test = wrf_e_ind(0) - wrf_s_ind(0) + 1
   test2 = wrf_e_ind(1) - wrf_s_ind(1) + 1

   wrf_prec = new((/n_years-1,n_hr_season/),float)
;  k is loop counter
   k = 1
   do i = 0, n_years - 2
        temp = dimsizes(f_in[:]->PREC_ACC_NC(wrf_s_ind(i):wrf_e_ind(i),nm1,nm2))
        wrf_prec(i,0:temp-1) = f_in[:]->PREC_ACC_NC(wrf_s_ind(i):wrf_e_ind(i),nm1,nm2)
        delete(temp)
        print(k)
        k = k + 1
   end do

;  Make data weekly totals
;  Convention: Week 1 begins at 1st value
   hr_wk = 24 * 7
   wrf_prec_wk = ndtooned(wrf_prec)
   temp = ceil(dimsizes(wrf_prec_wk)/hr_wk)
   n_wks = toint(temp)
   delete(temp)
   n_days = n_wks * 7
   weeks = ispan(1,281,1)
   wrf_prec_wks = new(n_wks-1,float)
   m = 1
   do i = 0, n_wks - 2
       wrf_prec_wks(i) = sum(wrf_prec_wk(hr_wk*i : (hr_wk*m)-1))
       m = m + 1
   end do

   wrf_prec_days = new(n_days-1,float)
   v = 1
   do i = 0, n_days - 2
       wrf_prec_days(i) = sum(wrf_prec_wk(24*i : (24*v)-1))
       v = v + 1
   end do

;    Import SNOTEL data to calculate sums to compare
   dirtxt = "/wind01/zrieck/CloudSeeding/Snotel/"
   filtxt = "Trial_Lake-daily.csv"
   pthtxt = dirtxt+filtxt
   strs   = asciiread(pthtxt,-1,"string")

   delim = ",-"
   nfields = str_fields_count(strs(1), delim)
;  print(nfields)

   iStrt = 1
   strend= ind(strs.eq.",,,,,,")
   iLast = strend(0)

   yyyy = toint(str_get_field(strs(iStrt:iLast), 1, delim))
   mm   = toint(str_get_field(strs(iStrt:iLast), 2, delim))
   dd   = toint(str_get_field(strs(iStrt:iLast), 3, delim))
   hhmn = toint(str_get_field(strs(iStrt:iLast), 4, delim))
   hh   = hhmn/100
   mn   = hhmn-(hh*100)
   prc  = tofloat(str_get_field(strs(iStrt:iLast),5, delim))
   qstr = str_get_field(strs(iStrt:iLast),6, delim)

   sec  = conform(yyyy, 0, -1)
   units= "Days since 2000-11-01 00:00:0.0"
   SN_time = cd_inv_calendar(yyyy,mm,dd,hh,mn,sec, units, 0)
   SN_time!0 = "SN_time"

;  Create grid for dataset
   grid2 = ispan(0,14600,1)
   grid2@units = "Days since 1978-10-01 00:00:0.0"
   grid_t2 = cd_calendar(grid2,-3)
;  Set data to skip Oct- giving Nov 1- March 31 every year
   sn_s_ind = new(n_years-1,integer)
   sn_e_ind = new(n_years-1,integer)
   do i = 0, n_years - 2
        temp1 = years(i) + "110100"
        temp2 = years(i+1) + "033100"
        sn_s_ind(i) = ind(grid_t2.EQ.temp1)
        sn_e_ind(i) = ind(grid_t2.EQ.temp2)
        delete(temp1)
        delete(temp2)
   end do

;  Manually specifying days between 110100 and 033113 since Apr-Oct is skipped
;  Leap years in 2004, 2008, and 2012
   Nov_days = 30
   Dec_days = 31
   Jan_days = 31
   Feb_days = 28
   Mar_days = 31
   n_lyr = 3
   tot_days = (Nov_days + Dec_days + Jan_days + Feb_days + Mar_days) * 13 + n_lyr
   tot_days_season = Nov_days + Dec_days + Jan_days + Feb_days + Mar_days
   n_wks_season_k = ceil(tot_days_season/7)
   n_wks_season = toint(n_wks_season_k)
   lyr = (/0,0,0,1,0,0,0,1,0,0,0,1,0/)
   prec = new(tot_days-1,float)
   do i = 0, n_years - 2
        q = tot_days_season + lyr(i)
        prec(q*i:q*(i+1)-1) = prc(sn_s_ind(i) : sn_e_ind(i))
   end do

;  Change precip data to weekly
   prec_wk = new(n_wks-1,float)
   d_wk = 7
   r = 1
   do i = 0, n_wks - 3
       prec_wk(i) = sum(prec(d_wk*i:(d_wk*r)-1))
       r = r + 1
   end do

;  Organize both plots into arrays separated by season
   wrf_wk_season = new((/n_years-1,n_wks_season/),float)
   y = 1
   do i = 0, n_years - 2
        wrf_wk_season(i,0:n_wks_season-1) = wrf_prec_wks((n_wks_season * i):(n_wks_season * y) - 1)
        y = y + 1
        print(y)
   end do

   sn_wk_season = new((/n_years-1,n_wks_season/),float)
   z = 1
   do i = 0, n_years - 2
       sn_wk_season(i,0:n_wks_season-1) = prec_wk((n_wks_season * i):(n_wks_season * z) -1)
        z = z + 1
   end do

;  Create an array with averages of every week through a season
   wrf_wk_avg = new(n_wks_season,float)
   do i = 0, n_wks_season - 1
       wrf_wk_avg(i) = avg(wrf_wk_season(:,i))
   end do

   sn_wk_avg = new(n_wks_season,float)
   do i = 0, n_wks_season - 1
       sn_wk_avg(i) = avg(sn_wk_season(:,i))
   end do

   print(wrf_wk_avg)
   print(sn_wk_avg)

   weeks_avg = ispan(0,n_wks_season-1,1)

;  1 plot will utilize cumulative sums
   prec_wk_cum = cumsum(sn_wk_avg,2)
   wrf_prec_wks_cum = cumsum(wrf_wk_avg,2)

;  Create 2d arrays for both plots
   scatter_data = new((/2,dimsizes(prec_wk)/),float)
   scatter_data(0,:) = prec_wk
   scatter_data(1,:) = wrf_prec_wks

   line_data = new((/2,dimsizes(prec_wk_cum)/),float)
   line_data(0,:) = prec_wk_cum
   line_data(1,:) = wrf_prec_wks_cum

;  Create scatter plot
   wks = gsn_open_wks("pdf", "Brighton_Test")
   res = True
   res@tiMainString = "Brighton Weekly Scater"
   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,prec_wk_cum,wrf_prec_wks_cum,res)
   weeks_avg = ispan(0,n_wks_season-1,1)

;  Create line plot
   wks2 = gsn_open_wks("pdf", "Brighton_Weekly_Cumulative")
   res2 = True
   colors = (/"blue","red"/)
   res2@xyLineColors = colors
   res2@tiMainString = "Brighton Cumulative Weekly Comparrison"
   res2@tmXBValues = (/ispan(0,20,5)/)
   res2@tmXBLabels = (/"Week 1","Week 6","Week 11","Week 16","Week 21"/)
   res2@tiXAxisString = "Number of Weeks"
   res2@tiYAxisString = "Precip (mm)"
   plot2 = gsn_csm_xy(wks2,weeks_avg,line_data,res2)

;  Insert label bar
   labels = (/"WRF","SNOTEL"/)
   lbres                    = True
   lbres@vpWidthF           = 0.1
   lbres@vpHeightF          = 0.1
   lbres@lbBoxMajorExtentF  = 0.36
   lbres@lbFillColors       = (/"red","blue"/)
   lbres@lbMonoFillPattern  = True
   lbres@lbLabelFontHeightF = 0.015
   lbres@lbLabelJust        = "CenterLeft"
   lbres@lbPerimOn          = False
   lbres@lgPerimColor   = "white"
   gsn_labelbar_ndc(wks2,2,labels,0.52,0.1,lbres)

;  Calculate total number of precipitation days when some amount of precipitation fell for each season
   wrf_days = new(n_years-1,integer)
   w = 1
   do i = 0, n_years - 2
       wrf_days(i) = sum(where(wrf_prec_days((i*tot_days_season):(w*tot_days_season)-1).GT.0.5,1,0))
       w = w + 1
   end do
   sn_days = new(n_years-1,integer)
   do i = 0, n_years - 2
        sn_days(i) = sum(where(prc(sn_s_ind(i):sn_e_ind(i)).NE.0,1,0))
   end do

    sres2 = True
    sres2@gsnDraw = True
    sres2@gsnFrame = False
    sres2@vpWidthF = 0.7
    sres2@vpHeightF = 0.5
    sres2@vpXF = .15
    sres2@trXMinF = 0
    sres2@trXMaxF = 13.5
    sres2@trYMinF = 0
    sres2@trYMaxF = 100
    sres2@tiXAxisString = "Season"
    sres2@tiYAxisString = "Frequency"
    sres2@tiMainString = "Brighton Precip Days Per Season"
    sres2@gsnXYBarChart = True
    sres2@gsnXYBarChartBarWidth = 0.10
    sres2@tmXBMode          = "Explicit"
    sres2@tmXBLabelFontHeightF = 0.010
    sres2@tmXTLabelFontHeightF = 0.010
    sres2@tmYLLabelFontHeightF = 0.010
    sres2@tiMainFontHeightF = 0.02
    sres2@tiMainFont = "helvetica"
    sres2@tmXBValues = (/ispan(1,13,1)/)
    sres2@tmXBLabels = (/"2000-01","2001-02","2002-03","2003-04","2004-05","2005-06","2006-07","2007-08","2008-09","2009-10","2010-11","2011-12","2012-13"/)

    plots3 = new((/1,2/),graphic)
    wks3 = gsn_open_wks("pdf","Brighton_Prec_Days")
    sres2@gsnXYBarChartColors = "red"
    plots3(0,0) = gsn_csm_xy(wks3,fspan(.675,12.675,13),wrf_days,sres2)
    sres2@gsnXYBarChartColors = "blue"
    plots3(0,1) = gsn_csm_xy(wks3,fspan(.925,12.925,13),sn_days,sres2)
    overlay(plots3(0,0),plots3(0,1))

;  Insert label bar
   gsn_labelbar_ndc(wks3,2,labels,0.52,0.1,lbres)

;   Output variables to Netcdf
;    cd "/wind01/zrieck/CloudSeeding/precip/Compilations"
    system("/bin/rm -f Data2.nc")
    ncdf = addfile("BR_Data2.nc" ,"c")
    filedimdef(ncdf,"time",-1,True)

;   Variables
;    ncdf->time  = f_time
    ncdf->BR_wrf_prec_cum = wrf_prec_wks_cum
    ncdf->BR_Snotel_prec_cum = prec_wk_cum
    ncdf->lat = a
    ncdf->lon = b

exit 
end