;;   For data import for Feb 10-19 2013
begin

;  Import WRF data

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

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

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

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

;---use gridpoint data
;  set lat-lon coordinates to match Cedar Creek
   a  = 41.90
   b  = -106.20
   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)

   sTime = tostring(f_in[:]->Times)

   yyyy_wrf = toint(str_get_cols(sTime, 0,3))
   mm_wrf  = toint(str_get_cols(sTime, 5,6))
   dd_wrf   = toint(str_get_cols(sTime, 8,9))
   hh_wrf   = toint(str_get_cols(sTime, 10,11))

;  Choose only Febrary data
   iYr = ind(yyyy_wrf.EQ.2013)
   iFeb = ind(mm_wrf.EQ.2)
   iDays = ind(dd_wrf.EQ.10)
  
   wrf_ind1 = venn2_intersection(iFeb,iYr)
   wrf_ind2 = venn2_intersection(iDays,wrf_ind1)
   wrf_s = wrf_ind2(17)
   wrf_e = wrf_ind2(17) + 215
   wrf_slw_f = f_in[:]->totalLiq(wrf_s:wrf_e,nm1,nm2)

;  Import surface observation data
   dirtxt = "/wind01/zrieck/CloudSeeding/slw/Obs/"
   filtxt = "CedarCreek_Feb10_19_2013.csv"
   pthtxt = dirtxt+filtxt
   strs   = asciiread(pthtxt,-1,"string")

   delim = ","
   test = str_split_csv(strs, ",", 0)

   nfields = str_fields_count(strs(1), delim)

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

   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))
   hh = toint(str_get_field(strs(iStrt:iLast), 4, delim))
   mn   = toint(str_get_field(strs(iStrt:iLast), 5, delim))
   slw  = tofloat(str_get_field(strs(iStrt:iLast),7, delim))
   qstr = str_get_field(strs(iStrt:iLast),6, delim)

   sec  = conform(yyyy, 0, -1)
   units= "Hours since 2013-02-10 17:00:0.0"
   ob_time = cd_inv_calendar(yyyy,mm,dd,hh,mn,sec, units, 0)
   ob_time!0 = "ob_time"

;  Find max values of SLW for every hour, and use those values as hourly SLW value
;  Create array of indices for every hour
   n_hr = 24
   hr23 = ind(hh.EQ.23)
   
   hrs = new((/n_hr,243/),integer)

   do i = 0, n_hr - 2
        hrs(i,:) = ind(hh.EQ.i)
   end do
   hrs(23,0:193) = hr23
 
;  Create array with indices for every day
   day = ispan(10,19,1)
   n_days = dimsizes(day)
   dimms = new(n_days,integer)
   j = 10
   do i = 0, n_days - 1
        dimms(i) = dimsizes(ind(dd.EQ.j))
        j = j + 1
   end do

   days = new((/n_days,dimms(1)/),integer)
   k = 10
   do i = 0, n_days - 1
       days(i,0:dimms(i)-1) = ind(dd.EQ.k)
       k = k + 1
   end do

;  Get indices for every hour of each day and find the single maximum SLW over that time
;  Test section
   test5 = venn2_intersection(hrs(23,:),days(1,:))
 
;  217 hours exist over data set with SOME data in them
   tot_hr = 217
   dimms3 = new((/1,n_hr/),integer)

   do i = 0, n_hr - 1
       dimms3(0,i) = dimsizes(venn2_intersection(hrs(i,:),days(0,:)))
   end do

   hr_days = n_hr * n_days
   indices = new((/hr_days,27/),integer)
   slw_f = new((/1,tot_hr-1/),float)

;  Day 0- only 7 hours starting at 17
   z = 0
   do i = 17, 23 
      temp = dimsizes(venn2_intersection(hrs(i,:),days(0,:))) 
      indices(i,0:temp-1) = venn2_intersection(hrs(i,:),days(0,:))
      slw_f(0,z) = max(slw(indices(i,0:temp-1)))
      z = z + 1
   end do

;  Day 1
   m = 0
   ma = 7 + (24 * 0)
   do i = n_hr, (n_hr * 2) - 1
       temp = dimsizes(venn2_intersection(hrs(m,:),days(1,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(m,:),days(1,:))
       slw_f(0,ma) = max(slw(indices(i,0:temp-1)))
       m = m + 1
       ma = ma + 1
       delete(temp)
   end do

;  Day 2
   n = 0
   na = 7 + (24 * 1)
   do i = (n_hr * 2), (n_hr * 3) - 1
       temp = dimsizes(venn2_intersection(hrs(n,:),days(2,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(n,:),days(2,:))
       slw_f(0,na) = max(slw(indices(i,0:temp-1)))
       n = n + 1
       na = na + 1
       delete(temp)
   end do

;  Day 3
   p = 0
   pa = 7 + (24 * 2)
   do i = (n_hr * 3), (n_hr * 4) - 1
       temp = dimsizes(venn2_intersection(hrs(p,:),days(3,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(p,:),days(3,:))
       slw_f(0,pa) = max(slw(indices(i,0:temp-1)))
       p = p + 1
       pa = pa + 1
       delete(temp)
   end do

;  Day 4
   q = 0
   qa = 7 + (24 * 3)
   do i = (n_hr * 4), (n_hr * 5) - 1
       temp = dimsizes(venn2_intersection(hrs(q,:),days(4,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(q,:),days(4,:))
       slw_f(0,qa) = max(slw(indices(i,0:temp-1)))
       q = q + 1
       qa = qa + 1
       delete(temp)
   end do

;  Day 5
   r = 0
   ra = 7 + (24 * 4)
   do i = (n_hr * 5), (n_hr * 6) - 1
       temp = dimsizes(venn2_intersection(hrs(r,:),days(5,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(r,:),days(5,:))
       slw_f(0,ra) = max(slw(indices(i,0:temp-1)))
       r = r + 1
       ra = ra + 1
       delete(temp)
   end do

;  Day 6
   s = 0
   sa = 7 + (24 * 5)
   do i = (n_hr * 6), (n_hr * 7) - 1
       temp = dimsizes(venn2_intersection(hrs(s,:),days(6,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(s,:),days(6,:))
       slw_f(0,sa) = max(slw(indices(i,0:temp-1)))
       s = s + 1
       sa = sa + 1
       delete(temp)
   end do

;  Day 7
   t = 0
   ta = 7 + (24 * 6)
   do i = (n_hr * 7), (n_hr * 8) - 1
       temp = dimsizes(venn2_intersection(hrs(t,:),days(7,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(t,:),days(7,:))
       slw_f(0,ta) = max(slw(indices(i,0:temp-1)))
       t = t + 1
       ta = ta + 1
       delete(temp)
   end do

;  Day 8
   u = 0
   ua = 7 + (24 * 7)
   do i = (n_hr * 8), (n_hr * 9) - 1
       temp = dimsizes(venn2_intersection(hrs(u,:),days(8,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(u,:),days(8,:))
       slw_f(0,ua) = max(slw(indices(i,0:temp-1)))
       ua = ua + 1
       u = u + 1
       delete(temp)
   end do

;  Day 9- 17 hours
   v = 0
   va = 7 + (24 * 8)
   do i = (n_hr * 9), (n_hr * 9 + 17) - 1
       temp = dimsizes(venn2_intersection(hrs(v,:),days(9,:)))
       indices(i,0:temp-1) = venn2_intersection(hrs(v,:),days(9,:))
       slw_f(0,va) = max(slw(indices(i,0:temp-1)))
       v = v + 1
       va = va + 1
       delete(temp)
   end do

;  Create date and totals array
   grid = ispan(0,215,1)
   grid@units = "Hours since 2013-02-10 17:00:0.0"
   grid_t = cd_calendar(grid, -3)
 
   totals = new((/2,dimsizes(wrf_slw_f)/),float)
   totals(0,:) = wrf_slw_f
   totals(1,:) = slw_f(0,:)

;  Create plot
   wks = gsn_open_wks("pdf", "CedarCreek_SLW")
   res = True
   res@trYMinF = 0
   res@trYMaxF = 0.5
   colors = (/"red","blue"/)
   res@tmXBMode     = "Explicit"
   res@tmXBValues = (/ispan(7,199,24)/)
   res@tmXBLabels = (/"Feb 11","Feb 12","Feb 13","Feb 14","Feb 15","Feb 16","Feb 17","Feb 18","Feb19"/)
   res@xyLineColors   = colors
   res@tiYAxisString = "SLW (mm)"
   res@tiXAxisString = "Day"
   res@tiMainString = "Cedar Creek SLW Feb 2013"
   plot = gsn_csm_xy(wks,grid_t,totals,res)
exit
end