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

begin

; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/moduser/weather/WRF37/WPS/geo_em.d01.nc","r")

  info	= (/ " 1  Evergreen Needleleaf"		, \  ;n=0
  	     " 2  Evergreen Broadleaf"  	, \
	     " 3  Deciduous Needleleaf"		, \
             " 4  Deciduous Broadleaf" 		, \
             " 5  Mixed Forest"	      	      	, \
             " 6  Closed Shrublands"		, \
             " 7  Open Shrublands"	      	, \
             " 8  Woody Savannas"	      	, \
             " 9  Savannas"		      	, \
             "10  Grasslands"		      	, \
             "11  Permanent Wetlands"	      	, \
             "12  Croplands"		      	, \
             "13  Urban and Built-up"	      	, \
             "14  Cropland/Natural Veg"		, \
             "15  Snow and Ice"	   	      	, \
             "16  Barren/Sparsely Veg"		, \
             "17  Water"    	     		, \
             "18  Wooded Tundra"  		, \
	     "19  Mixed Tundra"			, \
	     "20  Barren Tundra"		/) ; n=19

; We generate plots, but what kind do we prefer?
;  type = "x11"
;  type = "pdf"
  type = "png"
;  type = "ps"
;  type = "ncgm"
 
  wks = gsn_open_wks(type,"plt_landuse_zoom_SG")
  gsn_define_colormap(wks,"MPL_jet")
  wks2 = gsn_open_wks(type,"plt_landmask_zoom_SG")
  gsn_define_colormap(wks2,"temp_diff_1lev")


; Set some Basic Plot options
    y_start = 127
    y_end   = 135
    x_start = 105
    x_end   = 115  

  res = True
  res@MainTitle = "REAL-TIME WRF"

  pltres = True	      						; Plot resources
  pltres@PanelPlot = True
  pltres@FramePlot = False
  pltres@gsnMaximize = True 
  pltres@gsnFrame = False
  pltres@gsnDraw = False
  pltres@ZoomIn = True
  pltres@Xstart = x_start
  pltres@Ystart = y_start
  pltres@Xend = x_end
  pltres@Yend = y_end

  mpres = True
  mpres@mpFillOn			= False
  mpres@mpOutlineOn			= False
  mpres@mpGridLineColor 		= "Grey"
  mpres@mpGridLineDashPattern 		= 16
  mpres@mpGridLineThicknessF 		= 1.0
  mpres@mpGridSpacingF			= 0.25
  mpres@ZoomIn 				= True        ; set up map info for zoomed area
  mpres@Xstart 				= x_start
  mpres@Ystart 				= y_start
  mpres@Xend 				= x_end
  mpres@Yend 				= y_end



; LANDUSE , SOIL & VEGETATION FRACTION

      opts = res
      opts@cnFillOn = True
      opts@cnRasterModeOn = True
      opts@gsnSpreadColors = True
      opts@gsnSpreadColorStart = 35
      opts@gsnSpreadColorEnd = -1
      opts@lbAutoManage = False
      opts@pmLabelBarDisplayMode = "Always"
      opts@lbOrientation = "Horizontal"
      opts@lbLabelStride = 1	;skip every other label
      opts@pmLabelBarWidthF = 0.85
      opts@pmLabelBarHeightF = 0.3
      

      if(isfilevar(a,"LU_INDEX"))
        land = wrf_user_getvar(a,"LU_INDEX",0)
	land_zoom = land(y_start:y_end,x_start:x_end)
;       opts@cnFillPalette = "radar_1"
	opts@ContourParameters = (/ 0., 20., 1. /)
	opts@lbLabelPosition = "Center"
	opts@lbLabelAlignment = "BoxCenters"
	opts@lbLabelFontHeightF = 0.012
	opts@lbLabelFontAspectF = 1.5
	numbers20 = (/"1","2","3","4","5","6","7","8","9","10","11",\
  		       "12","13","14","15","16","17","18","19","20"/)
        opts@lbLabelStrings = ispan(0,20,1)
        contour = wrf_contour(a,wks,land_zoom,opts)

	rtxt  = True
	rtxt@txJust = "CenterLeft"
	rtxt@txFontHeightF = 0.010
	n  = -1
      	   xx = 0.1
      	   do ncol = 0,4
      	      yy   = 0.15
	     do nrow = 0,3
	     	n 	= n+1
	   	gsn_text_ndc (wks,info(n),xx,yy,rtxt)
	   	yy   = yy - 3*rtxt@txFontHeightF
	     end do
	      xx = xx + 0.155
      	   end do
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
      end if
      
       map_res = True
       map_res@mpOutlineOn 	= True
       map_res@mpFillOn    	= False
       dir 			= "/home/moduser/weather/NCL_programs/asia_hindcast/SCS_shp/"
       filename 		= "SCSENC_EA200001to4_LandArea_v201206.shp"
       new_map = gsn_add_shapefile_polylines(wks,plot,dir+filename,map_res)
      
  draw(plot)
  frame(wks)


      if(isfilevar(a,"LANDMASK"))
        landmask = wrf_user_getvar(a,"LANDMASK",0)
	landmask_zoom = landmask(y_start:y_end,x_start:x_end)
	opts@ContourParameters = (/ 0., 1., 1. /)
        contour = wrf_contour(a,wks2,landmask_zoom,opts)
        plot = wrf_map_overlays(a,wks2,(/contour/),pltres,mpres)
      end if
      
       map_res = True
       map_res@mpOutlineOn 	= True
       map_res@mpFillOn    	= False
       dir 			= "/home/moduser/weather/NCL_programs/asia_hindcast/SCS_shp/"
       filename 		= "SCSENC_EA200001to4_LandArea_v201206.shp"
       new_map = gsn_add_shapefile_polylines(wks2,plot,dir+filename,map_res)
      
  draw(plot)
  frame(wks2) 

delete(opts)

end