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

begin
;---Read data
  a    = addfile("FRP.nc","r")
  data = short2flt(a->frpfire(0,:,:))   ; lat x lon

;---Region of interest
  minlat = 26
  maxlat = 38
  minlon = 72
  maxlon = 86

  data_dims = dimsizes(data)
  latitude2d = conform_dims(data_dims,data&latitude,0)
  longitude2d = conform_dims(data_dims,data&longitude,1)

  data_mask = where(latitude2d.ge.minlat.and.latitude2d.le.maxlat.and.\
                    longitude2d.ge.minlon.and.longitude2d.le.maxlon,data,data@_FillValue)
  copy_VarMeta(data,data_mask)

;---Debug prints
  printMinMax(data,0)
  printMinMax(data_mask,0)

;---Start the graphics
 wks = gsn_open_wks("png","masked_box_plots1")

  res                   = True

  res@gsnDraw           = False    ; will panel later
  res@gsnFrame          = False

  res@gsnAddCyclic      = False    ; turn off longitude cyclic point

  res@cnFillOn          = True     ; turn on contour fill
  res@cnFillPalette     = "WhBlGrYeRe"
  res@cnLinesOn         = False    ; turn off contour lines
  res@cnLineLabelsOn    = False    ; turn off line labels
  res@lbLabelBarOn      = False    ; will add labelbar in panel

;  res@cnFillMode        = "RasterFill"
;  res@cnRasterSmoothingOn = True

;---Fix the contor levels for all plots
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/0.01,0.02,0.03,0.04,0.05/)
 
;---Zoom in on region of interest
  res@mpMinLatF         = minlat-2
  res@mpMaxLatF         = maxlat+2
  res@mpMinLonF         = minlon-2
  res@mpMaxLonF         = maxlon+2

  res@mpFillOn          = False     ; will add shapefile outlines later
  res@mpOutlineOn       = False

  res@tiMainString      = "original data"
  res@gsnCenterString   = data@long_name
  res@gsnLeftString     = ""
  res@gsnRightString    = ""

;---Plot of original data before masking
  plot_orig = gsn_csm_contour_map(wks,data,res)

;---Plot of data after masking
  res@tiMainString      = "masked data"
  plot_mask     = gsn_csm_contour_map(wks,data_mask,res)

  lnres = True
  lnres@gsLineThicknessF = 2.0
  lnres@gsLineDashPattern = 2
  dum1 = gsn_add_polyline(wks,plot_orig,(/minlon,minlon,maxlon,maxlon,minlon/),\
                                        (/minlat,maxlat,maxlat,minlat,minlat/),lnres)
  dum2 = gsn_add_polyline(wks,plot_mask,(/minlon,minlon,maxlon,maxlon,minlon/),\
                                        (/minlat,maxlat,maxlat,minlat,minlat/),lnres)


  gres                 = True
  gres@gsLineColor     = "red"
  gres@gsLineThicknessF  = "8"

   lon = (/77,78,78.0,78.5,79,79.0,80.0,80,81.5,81.5,82.5,82.5,83.5,83.5,84,84.0,83.0,83.0,81.0,81.0,80.0,80.0,78.0,78.0,77.0,77.0/)
  lat = (/30,30,29.5,29.5,29,28.5,28.5,28,28.5,27.5,27.5,26.5,26.5,26.0,26,28.5,28.5,29.5,29.5,30.0,30.0,30.5,30.5,31.0,31.0,30.0/)

  dum = gsn_add_polyline(wks,plot_mask,lon,lat,gres)



;---Panel original data and masked data
  pres = True
  pres@gsnMaximize = True
  pres@gsnPanelLabelBar = True
  gsn_panel(wks,(/plot_orig,plot_mask/),(/1,2/),pres)

;---Add filled dots at grid locations and panel
  mkres = True
  mkres@gsMarkerSizeF = 3
  mkres@gsMarkerIndex = 16
  mkres@gsnCoordsAttach = True
;  gsn_coordinates(wks,plot_orig,data,mkres)
;  gsn_coordinates(wks,plot_mask,data_mask,mkres)
;  gsn_panel(wks,(/plot_orig,plot_mask/),(/1,2/),pres)
end