; is a basic NCL template for creating contours over maps using
; WRF data, and attaching shapefile outlines. This template plots one 
; timestep of HGT from a WRF file, and attaches shapefile outlines
; from the "USA_adm2.shp" and "CAN_adm2.shp" shapefiles, downloaded
; from http://www.gadm.org/country/
;----------------------------------------------------------------------
; Note: as of NCL V6.2.0, you don't need these three load commands
;----------------------------------------------------------------------
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/wrf/WRFUserARW.ncl"

begin

  filename1 = "/DK/dk20/veiga/WRF-OUTPUT/ERILANE/wrfout-SC"
  filename2 = "/DK/dk20/veiga/WRF-OUTPUT/ERILANE/wrfout-CC"

  a        = addfile(filename1,"r")
  b        = addfile(filename2,"r")


  aa  = a->RAINNC(144,:,:)                        ; Selecting the last time-step
  bb  = a->RAINC(144,:,:)                         ; Selecting the last time-step
  cc  = aa+bb

  dd  = b->RAINNC(144,:,:)                        ; Selecting the last time-step
  ee  = b->RAINC(144,:,:)                         ; Selecting the last time-step
  ff  = dd+ee


  ; Computing anomaly

  anom = ff - cc


  wks = gsn_open_wks("pdf","teste")
  gsn_define_colormap(wks,"OceanLakeLandSnow")   ; Change color map

  res                   = True                   ; Use basic options for this field
  res@cnFillOn          = True                  ; Create a color fill plot

  res@ContourParameters = (/ -50., 50., 2/)    ; Contour intervals 

  contour = wrf_contour(a,wks,anom,res)  

  pltres            = True   ; Set plot options
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove contours

  mpres             = True   ; Set map options
  mpres@mpOutlineOn = False
  mpres@mpFillOn    = False

;---Create the contours over the WRF map (nothing will be drawn yet).
  plot = wrf_map_overlays(a,wks,contour,pltres,mpres) 

;---Attach the shapefile polylines using files read off gadm.org/country.
  munic_shp_name    = "../../SHAPES-BRASIL/AMAZONAS/MUNICIPIOS/municipios_AM.shp"
  bairros_manaus    = "../../SHAPES-BRASIL/AMAZONAS/MANAUS/BAIRROS_MANAUS.shp"
;  shp_rios        = "../../SHAPES-BRASIL/AMAZONAS/RIOS/rio_permanente.shp"
  shp_igor        = "../../SHAPES-BRASIL/AMAZONAS/Amazonia_legal_hidro_pol/AMAZ_LEG_HIDRO_POL_250.SHP"

;  canada_shp_name = "CAN_adm/CAN_adm2.shp"

  lnres                  = True
  lnres@gsLineColor      = "gray25"
  lnres@gsLineThicknessF = 0.5   
  munic_id = gsn_add_shapefile_polylines(wks,plot,munic_shp_name,lnres)
  man_id   = gsn_add_shapefile_polylines(wks,plot,bairros_manaus,lnres)
;  rio_id = gsn_add_shapefile_polylines(wks,plot,shp_rios,lnres) 
  igo_id = gsn_add_shapefile_polylines(wks,plot,shp_igor,lnres)

  draw(plot)       ; This will draw the map and the shapefile outlines.
  frame(wks)       ; Advance the frame

end