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

begin
;---Read data
  a    = addfile("/home/kunal/rain.nc","r")
  
  data = a->rfc(0,:,:)
  
 ; data2 = a->SWTNTCLRCLN(0,:,:)
  
; data = (data1-data2)

;---Region of interest
  minlat = 5
  maxlat = 38
  minlon = 60
  maxlon = 98

;---These are the two shapefiles we'll use for masking.
  shp_filename_ind = "IND_adm/india_state.shp"


;---Mask the data based on Nepal outline in shapefile
  opt             = True
  opt@minlat      = minlat    ; this makes the masking go faster
  opt@maxlat      = maxlat
  opt@minlon      = minlon
  opt@maxlon      = maxlon
  opt@debug       = True
 

;---Mask the data based on Northern India outlines in shapefile
 ; opt@shape_var   = "NAME_1"


  data_mask = shapefile_mask_data(data,shp_filename_ind,opt)
  
;---Combine the two sets of data back into "data_mask"
  
  


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

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

  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@cnFillMode        = "RasterFill"
  res@cnRasterSmoothingOn = True
  res@lbLabelBarOn      = True    ; will add labelbar in panel
  res@lbBoxLinesOn =    False
  res@lbLabelFontAspectF = 1.5
  res@cnFillPalette = "MPL_seismic"
  res@pmLabelBarWidthF            = 0.54
  res@pmLabelBarHeightF           = 0.090 ;adjuct label bar
  res@pmLabelBarOrthogonalPosF    = 0.10   ; adjuct position of the label bar 
 

;---Fix the contor levels for all plots
 ; res@cnLevelSelectionMode = "ManualLevels"  ; at SEASONAL TOA
 ; res@cnMinLevelValF = -0.8
 ; res@cnMaxLevelValF = 0.8
;  res@cnLevelSpacingF = 1


 
;---Zoom in on region of interest
  res@mpMinLatF         = minlat
  res@mpMaxLatF         = maxlat
  res@mpMinLonF         = minlon
  res@mpMaxLonF         = maxlon

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

;---Titles and tickmarks
   res@gsnLeftString     = "" 
;  res@gsnRightString    = ""
;  res@tiMainString      = ""
;  res@pmTitleZone       = 1
  res@pmTickMarkDisplayMode = "Always"   ; better tickmarks
  res@pmTickMarkDisplayMode  = "Always"
  res@mpGridSpacingF         =  1
  res@mpGridLineThicknessF   = 4.0
 

  
;  res@tmXBLabelAngleF = 90  
;  res@tmXBLabelFontAspectF  = 1
  res@tmXBLabelFontHeightF  = 0.018
;  res@tmYLLabelFontAspectF  = 1
  res@tmYLLabelFontHeightF  = 0.018
  res@tmXBLabelFontQuality  = "High"
  res@tmBorderThicknessF    = 6

  
 


;---Plot of data after masking
  res@tiMainString      = ""
  res@tiMainFontHeightF = 0.020 

  plot_mask     = gsn_csm_contour_map(wks,data_mask,res)
  
;---Add shapefile outlines to both plots
  lnres = True
  lnres@gsLineThicknessF = 4.0

   id1 = gsn_add_shapefile_polylines(wks,plot_mask,shp_filename_ind,lnres)
   



;---Panel original data and masked data
  pres = True
  pres@gsnMaximize = True
;  pres@gsnPanelLabelBar = True
  pres@cnFillPalette              = "GMT_seis"
  


  gsn_panel(wks,(/plot_mask/),(/1,1/),pres)


end