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

;----------------------------------------------------------------------
; Iran Provinces shapefile
undef("add_outlines_to_map")
procedure add_outlines_to_map(wks,plot)
local fnames, colors, lnres, n
begin
;
; These two shapefiles contain waterways and administrative 
; geographic info. We will add them to the existing map using
; polylines.
;
;  fnames = (/"TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3","../2shapefile/iran_province"/) + ".shp"
  fnames = (/"./shp_ostan/iran_province"/) + ".shp"

;  colors = (/"LightBlue", "Tan"/)    ; water, administrative
  colors = (/"Tan"/)    ; water, administrative

  lnres                  = True      ; resources for polylines
  lnres@gsLineThicknessF = 2.0       ; 2x as thick

; Ehsan: Following loop is unnecessary
  do n=0,dimsizes(fnames)-1
    dumstr = unique_string("primitive")
    lnres@gsLineColor = colors(n)
    plot@$dumstr$ = gsn_add_shapefile_polylines(wks, plot, fnames(n), lnres)
  end do
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

begin
;Part1: Reading CSV files(rbscReport_edited5_NW_R24h_DelNul.csv, rbscReport_edited5_NW_R24h_Meteomanz.csv)
  print("Part1: Reading .csv files (rbscReport_edited5_NW_R24h_DelNul.csv)") ;(rbscReport_edited5_NW_R24h_Meteomanz.csv)
  fname_NW      = "rbscReport_edited5_NW_R24h_DelNul.csv"
  lines_NW      = asciiread(fname_NW,-1,"string")
  delim = ","

;---Read fields 1-5 from fname_Report
  station_id_NW = tointeger(str_get_field(lines_NW(1:),1,delim))
  icao_code_NW  =           str_get_field(lines_NW(1:),2,delim)
  st_name_NW    =           str_get_field(lines_NW(1:),3,delim)
  lat_NW        = tofloat(str_get_field(lines_NW(1:),4,delim))
  lon_NW        = tofloat(str_get_field(lines_NW(1:),5,delim))
  Vdate_NW      =         str_get_field(lines_NW(1:),7,delim)
  R24h_NW       = tofloat(str_get_field(lines_NW(1:),8,delim))
  print(station_id_NW+"  "+icao_code_NW+"  "+st_name_NW+"  "+lat_NW+"  "+lon_NW)
  nSta          = dimsizes(station_id_NW)

;---Print the information
  print("End Part1 Reading .csv files (rbscReport_edited5_NW_R24h_DelNul.csv)")

;-- define the workstation (plot type and name)
wks = gsn_open_wks("png","plot_unstructured_grid_camse")
;-- set resources

j=0
latn=new(22,float)
lonn=new(22,float)
R24hn=new(22,float)
do i=0, nSta-1
  if (Vdate_NW(i) .eq. "20170414") then
    latn(j)=lat_NW(i)
	lonn(j)=lon_NW(i)
	R24hn(j)=R24h_NW(i)
	j=j+1
  end if
end do

res = True
res@gsnMaximize = True
res@cnFillOn = True ;-- turn on contour fill
;res@cnFillPalette = "BlueWhiteOrangeRed" ;-- choose color map
res@tiMainString = "NCL Doc Example: Unstructured grid (CAM-SE)" ;-- title
res@tiMainFontHeightF = 0.02
;---Lat/lon arrays of curvilinear grid for overlaying on map
res@sfXArray = lonn;lon1d
res@sfYArray = latn;lat1d

       res@gsnDraw              = False    ; don't draw
       res@gsnFrame             = False    ; don't advance frame
       res@mpDataSetName         = "Earth..4"   ; 1 through 4
	   res@mpOutlineBoundarySets = "National"
       res@mpFillOn             = False
       res@mpMinLatF             = 33.3;min(prc&lat)-1;min(prc&latitude)
       res@mpMaxLatF             = 41.3;max(prc&lat)+1;max(prc&latitude)
       res@mpMinLonF             = 40.0;min(prc&lon)-1;min(prc&longitude)
       res@mpMaxLonF             = 52.5;max(prc&lon)+1;max(prc&longitude)
       res@mpCenterLonF          = 0.5*(41.0 + 51.5)
       res@cnFillPalette    = BlAqGrYeOrReVi200_modified;(::-1)      ;reverse the color map
       res@gsnAddCyclic         = False    ; data is regional
       res@cnLinesOn            = False       ; Turn off contour lines
       res@cnLineLabelsOn       = False       ; Turn off contour lines
       ;res@cnFillMode           = "CellFill"  ; The CellFill method for non-rectangular ...
	   res@cnFillMode           = "AreaFill"
       ;res@cnFillMode           = "RasterFill"         ; Raster Mode
	   ;res@cnRasterSmoothingOn  = True
       res@lbLabelBarOn         = True
	   res@lbLabelAutoStride    = True
       res@lbBoxEndCapStyle     = "TriangleBothEnds"
       res@lbLabelFontHeightF   = 0.0150              ; change font size
       res@mpShapeMode          = "FreeAspect"
       res@vpWidthF             = 0.8
       res@vpHeightF            = 0.4


;-- draw the contour map
plot = gsn_csm_contour_map(wks,R24hn,res)
add_outlines_to_map(wks,plot)

       resP = True
       resP@gsnMaximize         = True                ; make ps/eps/pdf large [no effect x11]
	   ;resP@gsnPanelFigureStrings = figure_strings
	   ;resP@gsnPanelFigureStringsFontHeightF = 0.01
	   resP@wkPSResolution = 3000
       ;resP@gsnPaperOrientation = "Portrait"          ; force portrait
       ;*resP@gsnPanelLabelBar    = True                ; add common colorbar
       ;*resP@lbLabelAutoStride   = True
       ;*resP@lbLabelFontHeightF  = 0.0150              ; change font size
       ;resP@txString        = "GPM Level 3 IMERGV04 Late Daily 0.1 x 0.1 Degree Precipitation";prc@long_name;"APHRODITE: RU: "+ymdPlot
       gsn_panel(wks,plot,(/1,1/),resP)               ; now draw as one plot



end