begin

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"

;   Import files
    file_path =  "/cloud01/Model/WRF/NCAR_13yr/2d_ctrl/precip/"
    fils = systemfunc("ls /cloud01/Model/WRF/NCAR_13yr/2d_ctrl/precip/")
    WRF_prefix = "wrf2d_d01_CTRL_PREC_ACC_NC_"

;   Specify winter dates so we only use winter precip files
    date1 = (/200010, 200101, 200110, 200201, 200210, 200301, 200310, 200401, 200410, 200501, 200510, 200601, 200610, 200701, 200710, 200801, 200810, 200901, 200910, 201001, 201010, 201101, 201110, 201201, 201210, 201301/)
    date2 = (/200012, 200103, 200112, 200203, 200212, 200303, 200312, 200403, 200412, 200503, 200512, 200603, 200612, 200703, 200712, 200803, 200812, 200903, 200912, 201003, 201012, 201103, 201112, 201203, 201212, 201303/)

;   Run through do loop that pulls in all the files we need for winter   
    n = dimsizes (date1)
    do i = 0,25,1
          WRF00 = addfile(file_path + WRF_prefix + date1(i) + "-" + date2(i) + ".nc", "r")
    end do 
    
    WRF_prec_x = WRF00->PREC_ACC_NC(0,:,:)
    ;Precip is output as mm/min, need to correct to more usable value
    WRF_prec = WRF_prec_x * 60            
    WRF_time = WRF00->Time
    n_time = dimsizes(WRF_time)
    n_prec = dimsizes(WRF_prec)
    grid = ispan(0,n_time,1)
    grid@units = "Days since 2000-11-01 00:00:0.0"
    grid_t = cd_calendar(grid,-3)
    s_date = 2000110100
    e_date = 2001033100
    s_time2 = ind(grid_t.EQ.s_date)
    e_time2 = ind(grid_t.EQ.e_date)

    s = day_of_year(2000,11,01)
    s_t = 365 - s
    e = day_of_year(2001,03,31)
    e_t = s_t + e

;   set lat-lon coordinates to match Trial Lake
    s_lat = 40.68
    e_lat = 42
    s_lon = -110.95
    e_lon = -100

    lat = WRF00->XLAT
    lon = WRF00->XLONG

;   2-D data
    WRF_prec@lat2d = lat
    WRF_prec@lon2d = lon

;   get x-y indices for lat-lon data, convert back to verify accuracy
    nm = getind_latlon2d(lat,lon,s_lat,s_lon)

;---use gridpoint data
    WRF_prec_gridpoint_x = WRF00->PREC_ACC_NC(s_time2:e_time2,nm(0,0),nm(0,1)) ; to get single grid point
    WRF_prec_gridpoint = WRF_prec_gridpoint_x * 1440
    time_f = grid_t(s_time2:e_time2)

    date00 = (/2000110100, 2001010100, 2001110100, 2002010100, 2002110100, 2003010100, 2003110100, 2004010100, 2004110100, 2005010100, 2005110100, 2006010100, 2006110100, 2007010100, 2007110100, 2008010100, 2008110100, 2009010100, 2009110100, 2010010100, 2010110100, 2011010100, 2011110100, 2012010100, 2012110100, 2013010100/)
    date002 = (/2000123100, 2001033100, 2001123100, 2002033100, 2002123100, 2003033100, 2003123100, 2004033100, 2004123100, 2005033100, 2005123100, 2006033100, 2006123100, 2007033100, 2007123100, 2008033100, 2008123100, 2009033100, 2009123100, 2010033100, 2010123100, 2011033100, 2011123100, 2012033100, 2012123100, 2013033100/)

;Create do loop to get precip values for multiple dates
   ; do i = 0,25
;        s_date00 = date00(i)
;        e_date00 = date002(i)
;        s_time00 = ind(grid_t.EQ.s_date00)
;        e_time00 = ind(grid_t.EQ.e_date00)
       ; WRF_prec_gridpoint00 = WRF00->PREC_ACC_NC(s_time00:e_time00,nm(0,0),nm(0,1))
   ; end do

 n_prec = dimsizes(WRF_prec)
    WRF_prec00 = sum(WRF_prec_gridpoint)
   print(WRF_prec00 + "mm" )
;   print(s_time00)

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

;SN_s_date = 0    ; units in time since start date
;SN_e_date = 149
;SN_s_time = (SN_time.EQ.SN_s_date)
;SN_e_time = closest_val(SN_time,SN_e_date)

SN_prc = prc(8068:8218)
SN_prec = sum(SN_prc)
;print(SN_s_time)
;print(SN_e_time)

;    do i = 0,25
;        SN_s_date00 = date00(i)
;        SN_e_date00 = date002(i)
;        SN_s_time00 = ind(grid_t.EQ.s_date00)
;        SN_e_time00 = ind(grid_t.EQ.e_date00)
;    end do

print(SN_prec)

 res = True
    res@tiMainString = "Trial Lake 40.68 N 110.95W elev 9945 ft" 
    sums = (/WRF_prec00,SN_prec/)
    wks = gsn_open_wks("pdf", "histo")
    plot1 = gsn_histogram(wks,sums,res)
   
   test = ind(grid_t.EQ.date00(2))
;   print(test)
    print(WRF00(0))
end