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/csm/contributed.ncl"

;----------------------------------------------------------------------
; Function to attach a labelbar to an existing map, given the colors
; and the levels. In this case, there better be one fewer labels
; than colors, since we're only labeling the interior edges.
;----------------------------------------------------------------------
function attach_labelbar(wks,map,colors,labels)
local nboxes, lbid, lbres, amres
begin
  nboxes = dimsizes(colors)

  lbres                    = True             ; labelbar resources
  lbres@vpWidthF           = 0.08             ; width of labelbar
  lbres@vpHeightF          = 0.60             ; height of labelbar
  lbres@lbPerimOn          = False            ; turn off perimeter
  lbres@lbOrientation      = "Vertical"       ; the default
  lbres@lbLabelAlignment   = "InteriorEdges"  
  lbres@lbFillColors       = colors           ; list of colors to use
  lbres@lbMonoFillPattern  = True             ; fill them all solid
  lbres@lbLabelFontHeightF = 0.010            ; label font height

  lbid = gsn_create_labelbar(wks,nboxes,labels,lbres)
;
; Create some annotation resources indicating how we want
; to attach the labelbar to the plot. Here, we are using
; the center-left of the labelbar as the point which we
; are going to position it, and then we use
; amOrthogonal/ParallelPosF to move it right and up.
;                                                                 
  amres                  = True
  amres@amJust           = "CenterLeft"
  amres@amParallelPosF   =  0.39          ; Move to the right
  amres@amOrthogonalPosF = -0.10          ; Move up

  annoid = gsn_add_annotation(map,lbid,amres)

  return(annoid)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------

begin  

  raw = addfile("1B11.20100831.72869.7.HDF","r")
;  a = addfile("../HEDAS/earl-2010/EARL07L.2010083012/ATMOS/wrfout_d02_2010-08-31_05:00:00.nc","r")
  lat = raw->Latitude
  lon = raw->Longitude
  val_low_raw =  raw->lowResCh
  val_high_raw = raw->highResCh
   printMinMax(val_low_raw,0)
   printMinMax(val_high_raw,0)
  ghz = (/ "10.65", "10.65", "19.35", "19.35", "21.30", "37.00", "37.00", "85.00", "85.00" /)
  trmmpol = (/ "V", "H", "V", "H", "V", "V", "H", "V", "H" /)

  printVarSummary(raw)
  printVarSummary(lat)
  printVarSummary(lon)
  printVarSummary(val_low_raw)
  printVarSummary(val_high_raw)
  npixlo = dimsizes(val_low_raw(0,:,0))
  npixel = dimsizes(val_high_raw(0,:,0))
  nscan = dimsizes(val_high_raw(:,0,0))

  scanStart = 1050
  scanEnd = 1100

  print(npixlo)
  print(npixel)
  print(nscan)
  print(max(lat))
  print(min(lat))
  print(max(lon))
  print(min(lon))

  ; do scan = 1000,1200
  ;  print(str_concat( (/ sprinti("%i",scan), ":", sprinti("%0.2i/",toint(raw->Month(scan))), sprinti("%0.2i/",toint(raw->DayOfMonth(scan))), sprinti("%0.4i",toint(raw->Year(scan))), ", ", sprinti("%0.2i",toint(raw->Hour(scan))), ":", sprinti("%0.2i ",toint(raw->Minute(scan))), sprintf("%4.2f ", min(lat(scan,:))), "-", sprintf(" %4.2f",max(lat(scan,:))), ",", sprintf("%4.2f ", min(lon(scan,:))), "-", sprintf(" %4.2f", max(lon(scan,:))) /) ))
  ; end do

  ; switch from nscan x pixlo x 7 to nscan x 7 x pixlo
  val_low = val_low_raw(nscan|:,nchanlo|:,npixlo|:)
  val_high = val_high_raw(nscan|:,nchanhi|:,npixel|:)

  nlchannels = dimsizes(val_low(0,:,0))
  nhchannels = dimsizes(val_high(0,:,0))

  printVarSummary(lat) 
  printVarSummary(lon)

  opts = True
  opts@mpDataBaseVersion = "HighRes"
  opts@mpOutlineDrawOrder = "PostDraw"   ; draw continental outline last
  opts@mpGeophysicalLineColor = "Black"
  opts@mpGridLineColor = "Black"
  opts@mpLimbLineColor = "Black"
  opts@mpNationalLineColor = "Black"
  opts@mpPerimLineColor = "Black"
  opts@mpUSStateLineColor = "Black"
  opts@mpFillOn = False                  ; turn off map fill
  opts@mpOutlineOn = True
  opts@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
  opts@mpLimitMode="Corners"             ; choose range of map

  opts@gsnLeftString        = ""
  opts@gsnRightString        = ""

  opts@gsnMaximize = True  
  opts@gsnDraw     = False
  opts@gsnFrame    = False

;---The -0.5 and the +2 add a margin around the plot area.
  opts@mpLeftCornerLatF = min(lat(scanStart:scanEnd,0:npixel-1:2))-0.5
  opts@mpLeftCornerLonF = min(lon(scanStart:scanEnd,0:npixel-1:2))-0.5
  opts@mpRightCornerLatF = max(lat(scanStart:scanEnd,0:npixel-1:2))+2
  opts@mpRightCornerLonF = max(lon(scanStart:scanEnd,0:npixel-1:2))+2
  
  type = "png"

;---Resources for markers
  mkres               = True
  mkres@gsMarkerSizeF = 10.  ; Make marker larger

;---Levels for coloring markers by
  levels  = ispan(80,320,10)
  labels  = tostring(levels)
  nlevels = dimsizes(levels)
  colors  = (/17,24,32,39,46,54,61,68,76,83,90,98,105,112,119,127,134,\
              141,149,156,163,171,178,185,193,200/)

;---Arrays to hold markers
  dum_low = new(nlevels*nlchannels,graphic)
  dum_hgh = new(nlevels*nlchannels,graphic)

;---Lat/lon arrays for placing markers
  lon1d = ndtooned(lon(scanStart:scanEnd,0:npixel-1:2))
  lat1d = ndtooned(lat(scanStart:scanEnd,0:npixel-1:2))

  ncount = 0
  do lchannel = 0, nlchannels - 1
    chnnumber = sprinti("%0.4i", lchannel+1)
    chn = sprinti("%i", lchannel+1)
    gh  = ghz(lchannel)
    opts@tiMainString         = str_concat( (/"TB for TRMM ch ",chn," (",gh," GHz ", trmmpol(lchannel), ")" /) )
    wks = gsn_open_wks(type, str_concat((/ "CRTM_earl_72869-",chnnumber, "-raw" /) ))

    filled_square = NhlNewMarker(wks,"y",35,0.,0.,1.,1.,0.)
    mkres@gsMarkerIndex = filled_square    ; 16 will give you a filled dot

    gsn_define_colormap(wks,"BlAqGrYeOrRevi200") 
    i = NhlNewColor(wks,0.7,0.7,0.7)              ; add gray to colormap
    bt = ndtooned(val_low(scanStart:scanEnd,lchannel,:)/100. + 100.)

;---Create the map, but don't draw it.
    map = gsn_csm_map(wks, opts)
;
; Loop through the levels and collect the data that falls
; between or above/below the given level(s).
;
    do nl=0,nlevels
      if(nl.eq.0) then
        ii = ind(bt.le.levels(0))
      else if(nl.eq.nlevels) then
        ii = ind(bt.gt.levels(nlevels-1))
      else
        ii = ind(bt.gt.levels(nl-1).and.bt.le.levels(nl))
      end if
      end if
      if(.not.any(ismissing(ii))) then
;---Set the color for this group of markers
        mkres@gsMarkerColor = colors(nl)
;---Attach markers to map
        dum_low(ncount) = gsn_add_polymarker(wks,map,lon1d(ii),lat1d(ii), mkres)
        ncount = ncount + 1
      end if
      delete(ii)
    end do

;---Now attach a labelbar, draw map, and advance frame.
    lbid_low = attach_labelbar(wks,map,colors,labels)
     maximize_output(wks,False)
;    draw(map)
;    frame(wks)
    delete(wks)
  end do

  opts@mpLimitMode = "Corners"
  opts@mpLeftCornerLatF = min(lat(scanStart:scanEnd,0:npixel-1))-0.5
  opts@mpLeftCornerLonF = min(lon(scanStart:scanEnd,0:npixel-1))-0.5
  opts@mpRightCornerLatF = max(lat(scanStart:scanEnd,0:npixel-1))+2
  opts@mpRightCornerLonF = max(lon(scanStart:scanEnd,0:npixel-1))+2

  delete(bt)
  delete(lat1d)
  delete(lon1d)

  lon1d = ndtooned(lon(scanStart:scanEnd,0:npixel-1))
  lat1d = ndtooned(lat(scanStart:scanEnd,0:npixel-1))

  ncount = 0
  do hchannel = 0, nhchannels - 1
    channel = hchannel+nlchannels+1
    chnnumber = sprinti("%0.4i", channel)
    chn = sprinti("%i", channel)
    gh  = ghz(channel-1)
    opts@tiMainString         = str_concat( (/"TB for TRMM ch ",chn," (",gh," GHz ", trmmpol(channel-1), ")" /) )
    wks = gsn_open_wks(type, str_concat((/ "CRTM_earl_72869-",chnnumber,"-raw" /) ))
    gsn_define_colormap(wks,"BlAqGrYeOrRevi200") 
    i = NhlNewColor(wks,0.7,0.7,0.7)              ; add gray to colormap

    bt = ndtooned(val_high(scanStart:scanEnd,hchannel,:)/100. + 100.)
;    bt = ndtooned(val_high(scanStart:scanEnd,hchannel,:))
    print("val_high: hchannel = " + hchannel)
    printMinMax(bt,0)

    map = gsn_csm_map(wks, opts)

;---Add a filled square to the current marker table
    filled_square = NhlNewMarker(wks,"y",35,0.,0.,1.,1.,0.)
    mkres@gsMarkerIndex = filled_square    ; 16 will give you a filled dot
;
; Loop through the levels and collect the data that falls
; between or above/below the given level(s).
;
    do nl=0,nlevels
      if(nl.eq.0) then
        ii = ind(bt.le.levels(0))
      else if(nl.eq.nlevels) then
        ii = ind(bt.gt.levels(nlevels-1))
      else
        ii = ind(bt.gt.levels(nl-1).and.bt.le.levels(nl))
      end if
      end if
      if(.not.any(ismissing(ii))) then
;---Set the color for this group of markers
        mkres@gsMarkerColor = colors(nl)
;---Attach markers to map
        dum_hgh(ncount) = gsn_add_polymarker(wks,map,lon1d(ii),lat1d(ii), mkres)
        ncount = ncount+1
      end if
      delete(ii)
    end do
;---Now attach a labelbar, draw map, and advance frame.
    lbid_hgh = attach_labelbar(wks,map,colors,labels)
     maximize_output(wks,False)
;    draw(map)
;    frame(wks)
    delete(wks)
  end do
end