;;  Comparrison of WRF and SNOTEL for Pine Creek Station
;   38.88N, 112.25 W, elev = 8790 ft
;   Recycle WRF scrip using different gridpoint

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_in = addfiles(f_list, "r")
   ListSetType(f_in, "cat")

   f_time = f_in[:]->Times
;  printVarSummary(f_time)

   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")

;---use gridpoint data
;   set lat-lon coordinates to match Pine Creek
    a  = 38.88
    b  = -112.25
    lat2d = f_in[:]->XLAT
    lon2d = f_in[:]->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)

;   import and standardize time
    WRF_time = f_in[:]->Time
    n_time = dimsizes(WRF_time)

;   Creating grid of arbitrary length 5200 to exceed number fo days from 2000110100 to 2013033100
    grid = ispan(0,5200,1)
    grid@units = "Days since 2000-10-01 00:00:0.0"
    grid_t = cd_calendar(grid,-3)

;   Calculate number of days that are skipped. Create array that skips 0 indices on 1st pass and difference every other time
    doy_s = day_of_year(2001,11,01)
    doy_e = day_of_year(2001,03,31)
    diff = doy_s - doy_e
    diff_t = (/0, diff, diff, diff, diff, diff, diff, diff, diff, diff, diff, diff, diff/)

;   Loop through years 2000-2013 to calculate yearly winter sums from Nov 1 to Mar 31
    years = ispan(2000,2013,1)
    n_years = dimsizes(years)
    mon_length = (/30,31,31,28,31/)
    n_months = dimsizes(mon_length)

;   Leap years occur in 2004, 2008, and 2012 within data set
;   Create an array with 1s for leap years and 0s for not
    leap_yr = (/0,0,0,1,0,0,0,1,0,0,0,1,0/)
    n_lyr = dimsizes(leap_yr)

       Nov_hr = mon_length(0)*24
       Dec_hr = mon_length(1)*24
       Jan_hr = mon_length(2)*24
       Mar_hr = mon_length(4)*24

    Feb_hr = new(n_years-1, "integer")
    do i = 0, n_years-2
       Feb_hr(i) = (mon_length(3)+leap_yr(i))*24
    end do

;  Loop through each month, creating monthly totals for each year
;  Explanation of s_time eq: grabs index in days, so multiply by 24 to
;  convert to hours. Then subtract to account for days between Apr 1
;  and Oct 31 that were not counted. Finally, account for the potential
;  for a leap_yr adding 24 extra hours to the year
   Nov = new(n_years-1, "float")
   do i = 0, n_years-2
       s_date7  = years(i) + "110100"
       s_time7 = (ind(grid_t.EQ.s_date7) * 24)  - (diff_t(i) * 24 * i) + (leap_yr(i) * 24)
       e_time7 = s_time7 + Nov_hr-1
       Nov(i) = sum(f_in[:]->PREC_ACC_NC(s_time7:e_time7,nm1,nm2))
       print(Nov)
   end do

  Dec = new(n_years-1, "float")
   do i = 0, n_years-2
       s_date8  = years(i) + "120100"
       s_time8 = (ind(grid_t.EQ.s_date8) * 24)  - (diff_t(i) * 24 * i) + (leap_yr(i) * 24)
       e_time8 = s_time8 + Dec_hr
       Dec(i) = sum(f_in[:]->PREC_ACC_NC(s_time8:e_time8,nm1,nm2))
       print(Dec)
  end do

  Jan = new(n_years-1, "float")
   ct4 = 1
   do i = 0, n_years-2
       s_date9  = years(ct4) + "010100"
       s_time9 = (ind(grid_t.EQ.s_date9) * 24)  - (diff_t(i) * 24 * i) + (leap_yr(i) * 24)
       e_time9 = s_time9 + Jan_hr
       Jan(i)= sum(f_in[:]->PREC_ACC_NC(s_time9:e_time9,nm1,nm2))
       print(Jan)
       ct4 = ct4 + 1
   end do

   Feb = new(n_years-1, "float")
   ct5 = 1
   do i = 0, n_years-2
       s_date10  = years(ct5) + "020100"
       s_time10 = (ind(grid_t.EQ.s_date10) * 24)  - (diff_t(i) * 24 * i) + (leap_yr(i) * 24)
       e_time10 = s_time10 + Feb_hr(i)
       Feb(i) = sum(f_in[:]->PREC_ACC_NC(s_time10:e_time10,nm1,nm2))
       ct5 = ct5 + 1
       print(Feb)
   end do

   Mar = new(n_years-1, "float")
   ct6 = 1
   do i = 0, n_years-2
       s_date11  = years(ct6) + "030100"
       s_time11 = (ind(grid_t.EQ.s_date11) * 24)  - (diff_t(i) * 24 * i) + (leap_yr(i) * 24)
       e_time11 = s_time11 + Jan_hr
       Mar(i) = sum(f_in[:]->PREC_ACC_NC(s_time11:e_time11,nm1,nm2))
       print(Mar)
       ct6 = ct6 + 1
   end do

     WRF_tot = new(n_years-1, "float")
     do i = 0, n_years-2
        WRF_tot(i) = sum(Nov(i)+Dec(i)+Jan(i)+Feb(i)+Mar(i))
     end do

;   Generate array that has values for each monthly total
;   Convention: WRFxx where xx is the end year in March between 2000 and 2013
    WRF01 = (/Nov(0),Dec(0),Jan(0),Feb(0),Mar(0)/)
    WRF02 = (/Nov(1),Dec(1),Jan(1),Feb(1),Mar(1)/)
    WRF03 = (/Nov(2),Dec(2),Jan(2),Feb(2),Mar(2)/)
    WRF04 = (/Nov(3),Dec(3),Jan(3),Feb(3),Mar(3)/)
    WRF05 = (/Nov(4),Dec(4),Jan(4),Feb(4),Mar(4)/)
    WRF06 = (/Nov(5),Dec(5),Jan(5),Feb(5),Mar(5)/)
    WRF07 = (/Nov(6),Dec(6),Jan(6),Feb(6),Mar(6)/)
    WRF08 = (/Nov(7),Dec(7),Jan(7),Feb(7),Mar(7)/)
    WRF09 = (/Nov(8),Dec(8),Jan(8),Feb(8),Mar(8)/)
    WRF10 = (/Nov(9),Dec(9),Jan(9),Feb(9),Mar(9)/)
    WRF11 = (/Nov(10),Dec(10),Jan(10),Feb(10),Mar(10)/)
    WRF12 = (/Nov(11),Dec(11),Jan(11),Feb(11),Mar(11)/)
    WRF13 = (/Nov(12),Dec(12),Jan(12),Feb(12),Mar(12)/)
    Months = (/200011,200012,200101,200102,200103/)
    yfrac = yyyymm_to_yyyyfrac(Months,0)

;    Import SNOTEL data to calculate sums to compare
dirtxt = "/wind01/zrieck/CloudSeeding/Snotel/"
 filtxt = "Pine_Creek-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"

;   Data file begins on February 1, 1985
grid2 = ispan(0,14600,1)
    grid2@units = "Days since 1985-02-01 00:00:0.0"
    grid_t2 = cd_calendar(grid2,-3)

;   Calculate Snotel sums
    SN_tot = new(n_years-1, "float")
    do i = 0, n_years-2
        s_date2 = years(i) + "110100"
        e_date2 = years(i+1) + "033100"
        s_time4 = ind(grid_t2.EQ.s_date2)
        e_time4 = ind(grid_t2.EQ.e_date2)
        prec := prc(s_time4:e_time4)
        SN_tot(i) = sum(prec)
    end do

;Need monthly Snotel sums
    SN_Nov = new(n_years-1, "float")
    do i = 0, n_years-2
       s_date3 = years(i) + "110100"
       e_date3 = years(i) + "113000"
       s_time5 = ind(grid_t2.EQ.s_date3)
       e_time5 = ind(grid_t2.EQ.e_date3)
       Nov_prec := prc(s_time5:e_time5)
       SN_Nov(i) = sum(Nov_prec)
    end do

    SN_Dec = new(n_years-1, "float")
    do i = 0, n_years-2
       s_date3 = years(i) + "120100"
       e_date3 = years(i) + "123100"
       s_time5 = ind(grid_t2.EQ.s_date3)
       e_time5 = ind(grid_t2.EQ.e_date3)
       Dec_prec := prc(s_time5:e_time5)
       SN_Dec(i) = sum(Dec_prec)
    end do

    SN_Jan = new(n_years-1, "float")
    cnt1 = 1
    do i = 0, n_years-2
       s_date3 = years(cnt1) + "010100"
       e_date3 = years(cnt1) + "013100"
       s_time5 = ind(grid_t2.EQ.s_date3)
       e_time5 = ind(grid_t2.EQ.e_date3)
       Jan_prec := prc(s_time5:e_time5)
       SN_Jan(i) = sum(Jan_prec)
       cnt1 = cnt1 +1
    end do

;   To account for leap years, finds index for Mar 1 every year and subtarcts 1 (data is daily)
    SN_Feb = new(n_years-1, "float")
    cnt2 = 1
    do i = 0, n_years-2
       s_date3 = years(cnt2) + "020100"
       e_date3 = years(cnt2) + "030100"
       s_time5 = ind(grid_t2.EQ.s_date3)
       e_time5 = ind(grid_t2.EQ.e_date3) - 1
       Feb_prec := prc(s_time5:e_time5)
       SN_Feb(i) = sum(Feb_prec)
       cnt2 = cnt2 + 1
    end do

    SN_Mar = new(n_years-1, "float")
    cnt3 =1
    do i = 0, n_years-2
       s_date3 = years(cnt3) + "030100"
       e_date3 = years(cnt3) + "033100"
       s_time5 = ind(grid_t2.EQ.s_date3)
       e_time5 = ind(grid_t2.EQ.e_date3)
       Mar_prec := prc(s_time5:e_time5)
       SN_Mar(i) = sum(Mar_prec)
       cnt3 = cnt3 + 1
    end do

    SN01 = cumsum((/SN_Nov(0),SN_Dec(0),SN_Jan(0),SN_Feb(0),SN_Mar(0)/),0)
    SN02 = cumsum((/SN_Nov(1),SN_Dec(1),SN_Jan(1),SN_Feb(1),SN_Mar(1)/),0)
    SN03 = cumsum((/SN_Nov(2),SN_Dec(2),SN_Jan(2),SN_Feb(2),SN_Mar(2)/),0)
    SN04 = cumsum((/SN_Nov(3),SN_Dec(3),SN_Jan(3),SN_Feb(3),SN_Mar(3)/),0)
    SN05 = cumsum((/SN_Nov(4),SN_Dec(4),SN_Jan(4),SN_Feb(4),SN_Mar(4)/),0)
    SN06 = cumsum((/SN_Nov(5),SN_Dec(5),SN_Jan(5),SN_Feb(5),SN_Mar(5)/),0)
    SN07 = cumsum((/SN_Nov(6),SN_Dec(6),SN_Jan(6),SN_Feb(6),SN_Mar(6)/),0)
    SN08 = cumsum((/SN_Nov(7),SN_Dec(7),SN_Jan(7),SN_Feb(7),SN_Mar(7)/),0)
    SN09 = cumsum((/SN_Nov(8),SN_Dec(8),SN_Jan(8),SN_Feb(8),SN_Mar(8)/),0)
    SN10 = cumsum((/SN_Nov(9),SN_Dec(9),SN_Jan(9),SN_Feb(9),SN_Mar(9)/),0)
    SN11 = cumsum((/SN_Nov(10),SN_Dec(10),SN_Jan(10),SN_Feb(10),SN_Mar(10)/),0)
    SN12 = cumsum((/SN_Nov(11),SN_Dec(11),SN_Jan(11),SN_Feb(11),SN_Mar(11)/),0)
    SN13 = cumsum((/SN_Nov(12),SN_Dec(12),SN_Jan(12),SN_Feb(12),SN_Mar(12)/),0)

    sums1 = new((/2,dimsizes(WRF01)/),float)
    sums1(0,:) = WRF01
    sums1(1,:) = SN01

    sums2 = new((/2,dimsizes(WRF01)/),float)
    sums2(0,:) = WRF02
    sums2(1,:) = SN02

    sums3 = new((/2,dimsizes(WRF01)/),float)
    sums3(0,:) = WRF03
    sums3(1,:) = SN03

    sums4 = new((/2,dimsizes(WRF01)/),float)
    sums4(0,:) = WRF04
    sums4(1,:) = SN04

    sums5 = new((/2,dimsizes(WRF01)/),float)
    sums5(0,:) = WRF05
    sums5(1,:) = SN05

    sums6 = new((/2,dimsizes(WRF01)/),float)
    sums6(0,:) = WRF06
    sums6(1,:) = SN06

    sums7 = new((/2,dimsizes(WRF01)/),float)
    sums7(0,:) = WRF07
    sums7(1,:) = SN07

    sums8 = new((/2,dimsizes(WRF01)/),float)
    sums8(0,:) = WRF08
    sums8(1,:) = SN08

    sums9 = new((/2,dimsizes(WRF01)/),float)
    sums9(0,:) = WRF09
    sums9(1,:) = SN09

    sums10 = new((/2,dimsizes(WRF01)/),float)
    sums10(0,:) = WRF10
    sums10(1,:) = SN10

    sums11 = new((/2,dimsizes(WRF01)/),float)
    sums11(0,:) = WRF11
    sums11(1,:) = SN11

    sums12 = new((/2,dimsizes(WRF01)/),float)
    sums12(0,:) = WRF12
    sums12(1,:) = SN12

    sums13 = new((/2,dimsizes(WRF01)/),float)
    sums13(0,:) = WRF13
    sums13(1,:) = SN13

    wks = gsn_open_wks("pdf", "Pine_Creek")
    plots = new(12,graphic)
    res          = True
    res@gsnDraw  = False                          ; don't draw
    res@gsnFrame = False                          ; don't advance frame
    res@cnInfoLabelOn = False                     ; turn off cn info label
    colors = (/"blue","red"/)
    res@tmXBMode     = "Explicit"
    res@tmXBValues   = (/"Nov","Dec","Jan","Feb","Mar"/)
    res@tmXBLabels   = res@tmXBValues
    res@xyLineColors   = colors
    res@tiXAxis = "Year"
    res@tiYAxis = "LPE (mm)"

    plots(0) = gsn_csm_xy(wks,yfrac,sums1,res)
    plots(1) = gsn_csm_xy(wks,yfrac,sums2,res)
    plots(2) = gsn_csm_xy(wks,yfrac,sums3,res)
    plots(3) = gsn_csm_xy(wks,yfrac,sums4,res)
    plots(4) = gsn_csm_xy(wks,yfrac,sums5,res)
    plots(5) = gsn_csm_xy(wks,yfrac,sums6,res)
    plots(6) = gsn_csm_xy(wks,yfrac,sums7,res)
    plots(7) = gsn_csm_xy(wks,yfrac,sums8,res)
    plots(8) = gsn_csm_xy(wks,yfrac,sums9,res)
    plots(9) = gsn_csm_xy(wks,yfrac,sums10,res)
    plots(10) = gsn_csm_xy(wks,yfrac,sums11,res)
    plots(11) = gsn_csm_xy(wks,yfrac,sums12,res)

    resP                    = True                 ; modify the panel plot
    resP@gsnPanelMainString = "Pine Creek 38.88 N 112.25W elev 8790 ft"     ; set main title
    gsn_panel(wks,plots,(/3,4/),resP)               ; now draw as one plot
exit
end