; 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
  T_Krig = readAsciiTable("Mean_Temperature_day_night.txt",4 , "float",1)
  temp = T_Krig(:,2) 
  index_g = ind(temp.ne.-999.)
  lat = T_Krig(index_g,0) 
  lon = T_Krig(index_g,1) 
  HourlyT2 = T_Krig(index_g,2) 

 figurename = "wrf_d05_OKCmicronet_UHII_hourly_smallerD_29"
  wks = gsn_open_wks("eps",figurename)
  gsn_define_colormap(wks,"cosam") ; select

  res                             = True
  res@gsnFrame                    = False   ; So we can draw markers
  res@gsnMaximize                 = True

  res@gsnSpreadColors             = True
  res@pmTickMarkDisplayMode       = "Always"

  res@trGridType                  = "TriangularMesh"    ; The default if you
                                                        ; have 1D data
  res@cnLineLabelPlacementMode    = "Constant"
  res@cnLineDashSegLenF           = 0.3

;  res@cnLevelSelectionMode        = "ManualLevels"
;  res@cnMinLevelValF              = 15         	; 15.25
;  res@cnMaxLevelValF              = 50		; 49.75
;  res@cnLevelSpacingF             = 0.25

  res@cnFillOn                    = True
  res@cnLinesOn                   = True
  res@cnLineLabelsOn              = True
  res@cnLevelFlags		  = new(139,"string")
  res@cnLevelFlags(:)		  = "NoLine"
  res@cnLevelFlags(0::20)         = "LineAndLabel"
  res@tmYLLabelStride = 2
  res@tmXBLabelStride = 2

  res@lbBoxLinesOn                = False
  res@gsnRightString    = "UHII,~S~o~N~C"
;  res@gsnLeftString    = "Avg at 22-05 CST"
  res@gsnLeftString    = "Avg at 09-17 CST"

  res@sfXArray                    = lon
  res@sfYArray                    = lat

  res@mpProjection                = "LambertConformal"

  res@mpLimitMode                 = "Corners"

  res@mpFillOn                    = False
  res@mpOutlineDrawOrder          = "PostDraw"
  res@mpFillDrawOrder             = "PreDraw"
  res@mpOutlineBoundarySets       = "GeophysicalAndUSStates"
  res@tmYRMajorOutwardLengthF = 0
  res@tmYLMajorOutwardLengthF = 0
  res@tmXBMajorOutwardLengthF = 0
  res@tmXBMinorOutwardLengthF = 0
   res@tmYRLabelsOn = False
   res@tmXTLabelsOn = False

  res@mpUSStateLineColor          = "Gray10"
  res@mpUSStateLineDashPattern    = 2

; fmap = addfile("wrfinput_d05.nc","r")
; WRF_map_c(fmap, res, 0) 
; res@tfDoNDCOverlay       = True
;  north= 173; 101 ; to 230
;  west = 80; 120 ; 26
;  right = 158
;  south = 103;61
;  lat_LU=fmap->XLAT(0,:,:)
;  lon_LU=fmap->XLONG(0,:,:)
;  res@mpRightCornerLonF =lon_LU(north,right);     -96.2171
;  res@mpRightCornerLatF = lat_LU(north,right)
;  res@mpLeftCornerLonF =lon_LU(south,west);  -98.46469
;  res@mpLeftCornerLatF = lat_LU(south,west)

  res@mpLambertMeridianF =	-100.719
  res@mpLambertParallel2F =	34.539
  res@mpLambertParallel1F =	34.539
  res@mpCenterLatF =	35.51027
  res@mpCenterLonF =	-97.54489
  res@mpRightCornerLonF =	-97.3266
  res@mpRightCornerLatF =	35.64397
  res@mpLeftCornerLonF =	-97.76929
  res@mpLeftCornerLatF =	35.34035

  res@mpDataBaseVersion  = "Ncarg4_1"
  res@mpDataSetName      = "Earth..4"
  res@mpOutlineBoundarySets = "AllBoundaries"
 print(res)
  map = gsn_csm_contour_map(wks,HourlyT2,res)


  getvalues map@contour
    "cnLevels" : levels
    "cnFillColors" : colors
  end getvalues
  num_distinct_markers = dimsizes(levels)+1        ; number of distinct markers
  printVarSummary(HourlyT2)
  dimsobs = dimsizes(HourlyT2)
  print(num_distinct_markers)
  lat_new = new((/num_distinct_markers,dimsobs/),float,-999)
  lon_new = new((/num_distinct_markers,dimsobs/),float,-999)
  do j = 0, num_distinct_markers-1
    if (j.eq.0) then
      indexes = ind(HourlyT2.lt.levels(0))
    end if
    if (j.eq.num_distinct_markers-1) then
      indexes = ind(HourlyT2.ge.max(levels))
     print("color "+j+" indexes="+indexes)
    end if
    if (j.gt.0.and.j.lt.num_distinct_markers-1) then
      indexes = ind(HourlyT2.ge.levels(j-1).and.HourlyT2.lt.levels(j))
    end if

    if (.not.any(ismissing(indexes))) then
      npts_range = dimsizes(indexes)   ; # of points in this range.
      lat_new(j,0:npts_range-1) = lat(indexes)
      lon_new(j,0:npts_range-1) = lon(indexes)
     print("color "+j+"npts_range="+npts_range)
    end if
    delete(indexes)            ; Necessary b/c "indexes" may be a different
                               ; size next time.
  end do ; j = 0, num_distinct_markers-1
  gsres               = True
;  gsres@gsMarkerIndex = 4          ; circle at first
  gsres@gsMarkerIndex = 16          ; circle at first

  do j = 0, num_distinct_markers-1
    if (.not.ismissing(lat_new(j,0)))
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.016
      gsres@gsMarkerColor      = "black"
;      gsres@gsMarkerColor      = "white"
      gsn_polymarker(wks,map,lon_new(j,:),lat_new(j,:),gsres)
       gsres@gsMarkerIndex = 16          ; Use filled dots for markers.
      gsres@gsMarkerSizeF = 0.011
      gsres@gsMarkerColor      = colors(j)
      gsn_polymarker(wks,map,lon_new(j,:),lat_new(j,:),gsres)
    end if
  end do


  frame(wks)   ; Now advance the frame.
  system("convert -trim -density 300 -resize 1000x600 "+figurename+".eps /nsftor/xhu/public_html/OKCmicronet/Spatial/200901/"+figurename+".png")
  system("rm "+figurename+".eps")



end