;----------------------------------------------------------------------
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 = (/"/home/taghizadeh/shp_ostan/iran_province"/) + ".shp"

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

  dumstr = unique_string("primitive")
  lnres@gsLineColor = "Brown"
  plot@$dumstr$ = gsn_add_shapefile_polylines(wks, plot, fnames, lnres)

end

;----------------------------------------------------------------------
; Add markers for cities to given map plot.
;----------------------------------------------------------------------
undef("add_places_to_map")
procedure add_places_to_map(wks,plot)
local f, stids, stnames, sticaos, stlats, stlons, stelvs, stdims, xres, mkres
begin
;---Read shapefile with "places" information.
  f = addfile("/home/taghizadeh/shp_stations/stations1.shp","r")

  stids   = f->stid
  stnames = f->stnm
  sticaos = f->ICAOC
  stlons  = f->stlon
  stlats  = f->stlat
  stelvs  = f->stelv
  stdims  = dimsizes(stnames)

;---Set up text resources to label random areas of interest.
  txres               = True
  txres@txFontHeightF = 0.010
  txres@txJust        = "topcenter" ; "topright"

;---Set up marker resources to mark random areas of interest.
  mkres               = True
  mkres@gsMarkerIndex = 16
  mkres@gsMarkerColor = "brown"
  mkres@gsMarkerSizeF = 7.

  icao_of_interest = (/ "OITR","OITT","OITL","OICS","OITZ","OIGG", \
						"OICC","OIHH","OIIK","OICI","OICK","MAAA", \
						"ALKK","MASA","OIII","OIAW","OIFS","OIFM", \
						"OIIS","OING","OISY","OIBB","OISS","OIYY", \
						"OIMB","OIMM","OIMN","OIKB","OIKK","OIZH" /)


;---Looking for interested icao codes
  icao_of_interest@_FillValue = "missing"
  do i=0,stdims-1
    if(any(sticaos(i).eq.icao_of_interest)) then
      ii = ind(icao_of_interest.eq.sticaos(i))
      icao_of_interest(ii) = "missing"
      dumstr = unique_string("primitive")

;---Attach the marker to the map.
      plot@$dumstr$ = gsn_add_polymarker(wks, plot, stlons(i), stlats(i), mkres)
      dumstr = unique_string("primitive")
;

;---Attach the text string to the map. 
      plot@$dumstr$ = gsn_add_text(wks, plot, "  " + sticaos(i), \
                                    stlons(i), stlats(i), txres)
      delete(ii)
    end if
  end do
end

;----------------------------------------------------------------------
begin
;---Open WRF output file
  fi	= "wrfout_d02_2018-02-18_00:00:00"
  di	= "/home/taghizadeh/wrfmaps/operational/wrfouts/"
  f		= addfile(di + fi,"r")

  times			= wrf_user_getvar(f,"times",-1) ; get times in the file
  ntimes		= dimsizes(times)         ; number of times in the file
  print(times)
  years			= str_get_cols(times,0,3)
  months		= str_get_cols(times,5,6)
  days			= str_get_cols(times,8,9)
  dow			= day_of_week(stringtointeger(years),stringtointeger(months),stringtointeger(days))
  dowc			= new(ntimes,string)

  pny           = str_get_cols(times(0),0,3)
  pnm           = str_get_cols(times(0),5,6)
  pnd           = str_get_cols(times(0),8,9)
  pnh           = str_get_cols(times(0),11,12)

;---Extracting day of week
 do it=0,ntimes-1
   if (dow(it) .eq. 0) then
	dowc(it) = "Sun"
 else if (dow(it) .eq. 1) then
	dowc(it) = "Mon"
 else if (dow(it) .eq. 2) then
	dowc(it) = "Tue"
 else if (dow(it) .eq. 3) then
	dowc(it) = "Wed"
 else if (dow(it) .eq. 4) then
	dowc(it) = "Thu"
 else if (dow(it) .eq. 5) then
	dowc(it) = "Fri"
 else
	dowc(it) = "Sat"
   end if
   end if
   end if
   end if
   end if
   end if
 end do

;---Read u, v.
 do it=0,ntimes-1,3
   U10   = wrf_user_getvar(f,"U10",0) ; u at 10 m
;   printVarSummary(U10)
   print("Working on Time "+it)
   itc = sprinti("%0.3i", it)
 
   wks = gsn_open_wks("png","U10WRF"+pny+pnm+pnd+pnh+"_"+itc)
 
   u_res				= True
   
   u_res@gsnDraw		= False
   u_res@gsnFrame		= False
   u_res@gsnMaximize	= True
   u_res@gsnStringFont  = 25
   u_res@gsnStringFontHeightF = 0.01
   u_res@gsnRightString	= ""
   u_res@gsnLeftString	= "U at 10 m "+"("+U10@units+")" + \
							" Valid Time "+dowc(it)+" " + \
							times(it)
   
   u_res@lbBoxEndCapStyle            = "TriangleBothEnds"
   u_res@pmLabelBarHeightF            = 0.06        ; Make labelbar less thick
   u_res@lbLabelFontHeightF          = 0.014
   u_res@pmLabelBarOrthogonalPosF    = 0.01

;--These three lines necessary for plotting in native WRF map projection
   u_res	= wrf_map_resources(f,u_res)
   u_res@tfDoNDCOverlay = "NDCViewport" ; True  ; These two resources are required when overlaying
   u_res@gsnAddCyclic	= False ; plots in native WRF map projection.
   
   u_res@cnFillOn		= True
   u_res@mpOutlineOn	= True
   
;---Fix the settings that wrf_map_resources used for map outline color and thicknesses
   u_res@mpOutlineBoundarySets       = "AllBoundaries"     ; gives us more outlines
   u_res@mpDataSetName               = "Earth..4"          ; better resolution, use "Earth..2" if don't need the resolution
   u_res@mpGeophysicalLineColor      = "navy"
   u_res@mpNationalLineColor         = "navy"
   u_res@mpGeophysicalLineThicknessF = 3.0
   u_res@mpNationalLineThicknessF    = 3.0
;   u_res@mpLandFillColor		= "white"
;Following two lines specify Iran area with different background color   
   u_res@mpDataBaseVersion   = "MediumRes"
   u_res@mpMaskAreaSpecifiers = (/"Iran"/)
   
   ucontour = gsn_csm_contour_map(wks,U10,u_res)
   
   add_outlines_to_map(wks,ucontour)
   add_places_to_map(wks,ucontour)
   
   draw(ucontour)
   frame(wks)
   
   delete(u_res)
   
 end do

end