load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

; Read in our ascii station data
  nrow =  73
  ncol = 26
  data = asciiread("data/StationDataColorado_amount.txt",(/nrow,ncol/),"float") 

; Take lats and lons from the station data
  lon = data(:,0)
  lat = data(:,1)

; Make new arrays for the station data and model data
  prec = new((/24,73/),float)
  pmodp = new((/24,73/),float)

; Loop through the station data to write it to our new arrays.
  j = 0
  do i=2,25
    prec(j,:) = data(:,i)
    j = j+1
  end do

; Add our files we made in the first step.
  a = addfile("data/prec_diurnal_jja_rnty.nc","r")
  modp = a->PREC

; Create a new array to hold the difference amount between the station and model grid cell
  pdiff = new((/24,73/),float)

  llres = True
  llres@returnInt = True

; Loop through each station to find the corresponding grid cell in the model data and
;  write to new array
  do l=0,nrow-1
    locij = wrf_user_ll_to_ij(a,lon(l),lat(l),llres)
    locij = locij-1
    locX = locij(0)
    locY = locij(1)
    pmodp(:,l) = modp(:,locY,locX)
  end do

; Fix our time zone issue and calculate the difference
  pmodpn = new((/24,73/),float)
  pmodpn(0:16,:) = pmodp(7:23,:)
  pmodpn(17:23,:) = pmodp(0:6,:)
  pdiff = pmodpn-prec

;-----Options-----
  arr = (/-4.,-3.,-2.,-1.,0.0,1.,2.,3.,4./)
  colors = (/30,70,85,130,239,239,182,191,202,237/)
  labels = new(dimsizes(arr)+1,string)
  narr = dimsizes(arr)
  ncolors = dimsizes(colors)
  labels2 = (/"-4.0","-3.0","-2.0","-1.0","0.0","1.0","2.0","3.0","4.0"/)


  do m=0,23 ; Start a loop to go through each hour
;------------------------------ 
; Create X and Y arrays to hold the points for each range and initialize
; them to missing values.  We want to use num_distinct_markers
; different colors, so we need num_distinct_markers sets of X and
; Y points.
    num_distinct_markers = dimsizes(arr)+1        ; number of distinct markers
    lat_new = new((/num_distinct_markers,dimsizes(pdiff(m,:))/),float,-999)
    lon_new = new((/num_distinct_markers,dimsizes(pdiff(m,:))/),float,-999)

; Group the points according to which range they fall in. At the
; same time, create the label that we will use later in the legend.
    do i = 0, num_distinct_markers-1
      if (i.eq.0) then
        indexes = ind(pdiff(m,:).lt.arr(0))
        labels(i) = "x < " + arr(0)
      end if
      if (i.eq.num_distinct_markers-1) then
        indexes = ind(pdiff(m,:).ge.max(arr))
        labels(i) = "x >= " + max(arr)
      end if
      if (i.gt.0.and.i.lt.num_distinct_markers-1) then
        indexes = ind(pdiff(m,:).ge.arr(i-1).and.pdiff(m,:).lt.arr(i))
        labels(i) = arr(i-1) + " <= x < " + arr(i)
      end if

; Now that we have the set of indexes whose values fall within 
; the given range, take the corresponding lat/lon values and store
; them, so later we can color this set of markers with the appropriate
; color.
      if (.not.any(ismissing(indexes))) then
        npts_range = dimsizes(indexes)   ; # of points in this range.
        lat_new(i,0:npts_range-1) = lat(indexes)
        lon_new(i,0:npts_range-1) = lon(indexes)
      end if
      delete(indexes)            ; Necessary b/c "indexes" may be a different
                                 ; size next time.
    end do

;=========================================================================
; Begin plotting section
    wks = gsn_open_wks("X11","precip_diff_hr"+ (m+1) + "_rnty")
    gsn_define_colormap(wks,"rainbow+white+gray")
    cmap = gsn_retrieve_colormap(wks)

    labelbarcolors = new((/ncolors,3/),float)
    do n = 0, ncolors-1
      labelbarcolors(n,:) = cmap(colors(n),:)
    end do

;Set up some map resources.
    mpres              = True
    mpres@gsnFrame     = False            ; Don't advance the frame
    mpres@mpDataBaseVersion     = "Ncarg4_1"
    mpres@mpFillOn              = True
    mpres@mpLandFillColor       = "Background"
    mpres@mpPerimOn             = True
    mpres@mpPerimDrawOrder      = "PostDraw"
    mpres@mpUSStateLineThicknessF = 3.0
    mpres@mpOutlineOn            = True
    mpres@mpOutlineBoundarySets  = "AllBoundaries"

; Zoom in on Colordao
    mpres@mpLimitMode       = "LatLon"            ; choose range of map
    mpres@mpMinLatF  = 35.5
    mpres@mpMinLonF  = -110
    mpres@mpMaxLatF = 42.5
    mpres@mpMaxLonF = -101
    mpres@mpCenterLonF = -105

    mpres@tiMainFontHeightF = 0.03
    mpres@tiMainString = "Precip Diff Model-Obs Hour: " + (m+1)
    map = gsn_csm_map(wks,mpres)

;  Create logical variables to hold the marker and text resources.
;  These markers are different than the XY markers, because they are not
;  associated with an XY plot. You can put these markers on any plot.
    gsres               = True
    gsres@gsMarkerIndex = 16          ; Use filled dots for markers.

    txres               = True
    txres@txFontHeightF = 0.01

; Loop through each grouping of markers, and draw them one set at
; a time, assigning the proper color and size with gsn_marker.
    do k=0,ncolors-1
      if(.not.ismissing(lat_new(k,0))) then
        gsres@gsMarkerColor = colors(k)
        gsn_polymarker(wks,map,lon_new(k,:),lat_new(k,:),gsres)
      end if
    end do

;*************************************************
; add custom label bar
;************************************************* 
    lbres = True
    lbres@lbPerimOn            = False               ; no label bar box
    lbres@lbOrientation        = "Vertical"        ; orientation
    lbres@vpWidthF             = 0.1                 ; size
    lbres@vpHeightF            = 0.42
    lbres@lbLabelFontHeightF   = 0.012               ; label font height
    lbres@lbLabelAlignment     = "InteriorEdges"     ; where to label
    lbres@lbMonoFillPattern    = True                ; fill sold
    lbres@lbFillColors         = labelbarcolors      ; must be RGB triplets
    txres = True
    txres@txFontHeightF = 0.015

    gsn_labelbar_ndc (wks,narr+1,labels2,0.77,0.71,lbres)
    gsn_text_ndc(wks,"mm",0.81,0.27,txres)

    frame(wks)    ; Advance the frame.

  end do

end