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("dataforfigure2.nc","r")
  data = a->MYD08_D3_6_Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean(0,:,:)   ; lat x lon

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

  data_dims = dimsizes(data)
  lat2d = conform_dims(data_dims,data&lat,0)
  lon2d = conform_dims(data_dims,data&lon,1)

  data_mask = where(lat2d.ge.minlat.and.lat2d.le.maxlat.and.\
                    lon2d.ge.minlon.and.lon2d.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_plots")

  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.3,0.5,0.7,0.9,1.,2./)
 
;---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)

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