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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;-------creat palette function---------;
function precip_pal()
local cmap, ncolors
begin
  cmap = (/ (/86,86,86/),(/0,255,0/),(/255,0,0/),(/255,0,255/),(/255,255,255/) /)

;---Create RGBA array from RGB array so we can set a couple colors to transparent 
   ncolors          = dimsizes(cmap(:,0))
   cmap_rgba        = new((/ncolors,4/),float)
   cmap_rgba(:,0:2) = cmap/255.
   cmap_rgba(:,3)   = 1    ; set all colors to fully opaque
   cmap_rgba(0,3) = 0   ; set first and last colors to fully transparent
   cmap_rgba(4,3)   =0
  return(cmap_rgba)
end

;------------------
;main code
;------------------
begin
;;;;;;;;;;;1. set for the  plot;;;;;;;;;;;;;;;;
  disType   = "pdf"
  disName   ="TsEvaluation"
  wks       = gsn_open_wks(disType,disName)

  ;---set for the boundary of the plot
  LeftCornerLatF      = 13.5
  LeftCornerLonF      = 94.5
  RightCornerLatF     = 45.5
  RightCornerLonF     = 135.5
  
;;;;;;;;;2. read in the wrf data and get the hgt for ploting;;;;;;;;;;;;;;;;;;;;;;;;;;
  wrfDataDir     ="/public/home/GlobeDisRadi/wrfcase/HuaiHeJuly/RE140721_140727/wrfv3_in_out/"
  wrfFileName    ="wrfout_d01_2014-07-21_00_00_00"
  wrfFile        = wrfDataDir + wrfFileName
  wrf            =addfile(wrfFile,"r")
 
  lat2d      = wrf_user_getvar(wrf,"XLAT",0)           ; 	[south_north] x [west_east]
  lon2d      = wrf_user_getvar(wrf,"XLONG",0)
  nWrfLat2d  = dimsizes(lat2d)
  nWrfLat    = nWrfLat2d(0)  
  nWrfLon    = nWrfLat2d(1)  
  minLatWrf  = min(lat2d)      
  maxLatWrf  = max(lat2d)
  minLonWrf  = min(lon2d)
  maxLonWrf  = max(lon2d)
  
;*******3 read in trmm file for cordination ***********************
  trmmDataDir     ="/public/home/GlobeDisRadi/wrfcase/HuaiHeJuly/TRMM201407/TRMM/changed_names/"
  trmmFileName    ="3B42.20140721.00.7.nc"
  trmmFile        =trmmDataDir + trmmFileName
  trmm            =addfile(trmmFile,"r") 
  
  trmmlat = trmm->latitude({minLatWrf:maxLatWrf})
  trmmlon = trmm->longitude({minLonWrf:maxLonWrf}) 
  nTrmmLat    = dimsizes(trmmlat)
  nTrmmLon    = dimsizes(trmmlon)
  minTrmmLat  = min(trmmlat)
  maxTrmmLat  = max(trmmlat)
  minTrmmLon  = min(trmmlon)
  maxTrmmLon  = max(trmmlon)
;*******4. read in the ascii data of TS evaluation results*********
 TsFiles  = systemfunc("ls ./Ts_matrix*.txt")
 nTsFiles = dimsizes(TsFiles)
 Ts       = new((/nTsFiles,nTrmmLat,nTrmmLon/),"integer")
 do m =0, nTsFiles-1
   Ts(m,:,:) = asciiread(TsFiles(m),(/nTrmmLat,nTrmmLon/),"integer")
 end do 

;******4.1 exame the readin data **********
ts_line=new((/nTsFiles,nTrmmLat/),string)
do m =0, nTsFiles-1
  do i=0,nTrmmLat-1
    ts_line(m,i)=str_concat(sprinti ("%1i,", Ts(m,i,:)))
  end do
  exameName = "exameTs_" + tostring(m) + ".txt"
  asciiwrite(exameName,ts_line(m,:))
end do 

;******4.2 give the cordination************

lat            =fspan(minTrmmLat,maxTrmmLat,nTrmmLat)
lat!0          ="lat"
lat&lat        =lat
lat@units      = "degrees_north"
lat@long_name  = "Latitude"

lon=fspan(minTrmmLon,maxTrmmLon,nTrmmLon)
lon!0          ="lon"
lon&lon        =lon
lon@units      ="degrees_east"
lon@long_name  ="Longtitude"

Ts!1    ="lat"
Ts&lat  =lat
Ts!2    ="lon"
Ts&lon  =lon

;****** 5.plot the data *********************
;****** 5.1 set for the shared res ****************
res                 = True

res@gsnAddCyclic    = False
res@gsnMaximize     = False               
res@gsnDraw         = False
res@gsnFrame        = False
        
res                 = wrf_map_resources(wrf,res)
res@mpMinLatF       = LeftCornerLatF
res@mpMaxLatF       = RightCornerLatF
res@mpMinLonF       = LeftCornerLonF
res@mpMaxLonF       = RightCornerLonF

res@mpDataBaseVersion = "MediumRes"    
res@mpDataSetName = "Earth..4"         
res@mpOutlineSpecifiers = "China:states"
res@mpLandFillColor     = "white"
res@mpNationalLineColor         = "Black"
res@mpGeophysicalLineColor      = "Black" 
res@mpProvincialLineColor       = "Black"    
res@mpGeophysicalLineThicknessF = 1  ; changes the thickness of continental outlines.
res@mpProvincialLineThicknessF  = 1

;set the grid line
  res@mpGridAndLimbOn        = True     ; defalt is True
  res@mpGridLineColor        = "grey58"  ; the color is deeper, the number is bigger.
  res@mpGridLineDashPattern  = 11
  res@mpGridSpacingF         = 5
  res@mpGridLineThicknessF   = 1
  res@mpLimbLineColor        = "Black"
  res@mpPerimLineColor       = "Black"

  res@tmXBOn                = False 
  res@tmYLOn                = False
  res@tmBorderThicknessF    = 1.0

  res@lbLabelBarOn        = False 

res@cnLevelSelectionMode = "ExplicitLevels"
res@cnFillOn             = True
res@cnLinesOn            = False
res@cnLineLabelsOn       = False
res@cnFillMode           ="RasterFill"
res@cnLevels             =(/0.5,1.5,2.5,3.5/)
res@cnFillPalette       = precip_pal()

res@tfDoNDCOverlay  = False
plot_ts = new(nTsFiles,"graphic")
do m =0,nTsFiles-1
plot_ts(m) = gsn_csm_contour_map (wks,Ts(m,:,:),res)
end do 

panel_res = True
  panel_res@gsnPanelLabelBar = True
  panel_res@gsnPanelRowSpec  = True 

  gsn_panel(wks,plot_ts,(/3/),panel_res)

frame(wks)
;******************************************

end