load "./calculate_hourly_values.ncl"
;=========================================
;   Import surface observation data
   dirtxt = "/wind01/zrieck/CloudSeeding/slw/Obs/"
   dirtxt = "./"
   filtxt = "Cedar_Creek_Radiometer.csv"
   pthtxt = dirtxt+filtxt
   strs   = asciiread(pthtxt,-1,"string")
   nline  = dimsizes(strs)                  ; first line is a header

;============
; Eliminate 'e' at the end of a number. Not recognized by NCL. [Maybe Excel?]
; This is line 1503. Note the 'e'
; 2008,12,28,22,38,52,270.5,0,0,1.27150982867121e,5,,,
; 
   strs(1:) = str_sub_str(strs(1:),"e","")
;============

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

   nfields = str_fields_count(strs(1), delim)

   iStrt = 1
   iLast = nline-1  ; last line

   print("nfields="+nfields+"  iStrt="+iStrt+" iLast="+iLast)
   print("-----")
;;;iStrt = 1        ; debug   >= 1
   iLast = 1500     ; debug   <= nline

   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))
   qstr = str_get_field(strs(iStrt:iLast),6, delim)

;print(str_get_field(strs(iStrt:iLast), 8, delim)+"  "+ \
;      str_get_field(strs(iStrt:iLast), 9, delim)+"  "+ \
;      str_get_field(strs(iStrt:iLast),10, delim))
;print(str_get_field(strs(iStrt:iLast),10, delim))

   slw1 = tofloat(str_get_field(strs(iStrt:iLast), 8, delim))
   slw2 = tofloat(str_get_field(strs(iStrt:iLast), 9, delim))
   slw3 = tofloat(str_get_field(strs(iStrt:iLast),10, delim))

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

;  Only need 1 slw value. We will use the maximum of the 3 values for each time point
;  Creating this array

   n_slw = dimsizes(slw1)
   print("n_slw="+n_slw)
   print("-----")

   values = new((/3,n_slw/),float)
   values(0,:) = slw1
   values(1,:) = slw2
   values(2,:) = slw3
   slw_max = new(n_slw, float)
   slw_avg = new(n_slw, float)

   do i = 0, n_slw-1
       slw_max(i) = max(values(:,i))
       slw_avg(i) = avg(values(:,i))    ; original had ?typo?   max(values(:,i))
   end do

   slw_max!0 = "ob_time"
   slw_avg!0 = "ob_time" 
   slw_max&ob_time = ob_time
   slw_avg&ob_time = ob_time

   printVarSummary(slw_max)
   printMinMax(slw_max,0)
   print("------")
   printVarSummary(slw_avg)
   printMinMax(slw_avg,0)
   print("------")

                ; hourly avg of max values
   slwHourAvgMax = calculate_hourly_values(slw_max, "avg", 0, False)

   printVarSummary(slwHourAvgMax)
   printMinMax(slwHourAvgMax,0)
   print("------")

   ymdh = cd_calendar(slwHourAvgMax&ob_time, -3)
   print(ymdh+"  "+slwHourAvgMax)
exit

;  Data is by the minute, need to convert to hourly
;  Begins at 23 min past the hour on December 27 for hour 1900
   n_min = 60
   s_min = 23
   x_hours = floor((n_slw/60))
   n_hours = toint(x_hours)
   skip = n_min - s_min
   interval = skip + 60

;  Start on 1st full hour of data
   slw_start = slw_avg(skip)

;  Hourly conversion
   slw_hourly = new(n_hours-1,float)
   do i = 0, n_hours-2
       hr_strt = skip + i * n_min
       hr_end = interval + i * n_min
       slw_hourly(i) = avg(slw_avg(hr_strt:hr_end))
   end do