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/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/sti/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/sti/gmt2loc.ncl"

begin
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
  dir = "/home/users/operation_data/typhoon/AHW/11W/2014080412/wrfprd/"
  fs = systemfunc("ls " + dir + "wrfout_d01_*")
  a = addfile(fs(0) + ".nc","r")
  print("fs=" + fs(0))
  print("fs=" + fs)
  nfs= dimsizes(fs)
  print("nfs=" + nfs)

; We generate plots, but what kind do we prefer?
;  type = "x11"
;  type = "pdf"
  type = "ps"
;  type = "ncgm"
;  type@wkOrientation = "landscape"
  type@wkOrientation = "portrait"
  wks = gsn_open_wks(type,"rain24hd01zooml")
;  gsn_define_colormap(wks,"wh-bl-gr-ye-re")
; Set some Basic Plot options
  ARWres = True
  ARWres@MainTitle                   = "SHA-THRAPS Prediction"
  ARWres@Footer                      = False
  ;ARWres@lbTitleOn                   = False
  ARWres@pmLabelBarOrthogonalPosF   = -0.07   ; move bar up or close
        ARWres@gsnMaximize          = True
        ARWres@gsnPaperOrientation = "portrait" 
; add china map data dir
  dirdat = "/home/users/yangyh/wrfpp/CHINA_mp/NEW/"

    f = addfile(fs(0)+".nc","r")
    rain = wrf_user_getvar(f,"RAINNC",0)
    dims = dimsizes(rain)
    print("dims=  " + dims)
    raintot = new ((/nfs-1,dims(0),dims(1)/), typeof(rain))
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  gsn_define_colormap(wks,"wh-bl-gr-ye-re")

   do ifs=24,nfs-1,24
    f = addfile(fs(ifs)+".nc","r")
    f2 = addfile(fs(ifs-24)+".nc","r")
;    times(ifs) = wrf_user_list_times(f)
;      print("times(ifs)=" + times(ifs))
  ac     = new(nfs,integer)
  ac2     = new(nfs,integer)
;
;    mpres = True                   ; Create map background
    mpres = ARWres                   ; Create map background
    mpres@mpFillOn = False
    i = NhlNewColor(wks,0.7,0.7,0.7)              ; add gray to colormap
    mpres@mpLandFillColor       = "black"            ; set land to be gray
    mpres@mpOceanFillColor       = "white"
    mpres@mpGeophysicalLineColor      = "black"
    mpres@mpGeophysicalLineThicknessF = 1.8
    mpres@mpDataBaseVersion = "Ncarg4_1"
    map = wrf_map(wks,f,mpres)
     ;zoom 
     ip_latss = (/15.0/)
     ip_lonss = (/100.0/)
     ip_latse = (/30.0/)
     ip_lonse = (/133./)
     locs =  wrf_user_latlon_to_ij(f, ip_latss, ip_lonss)
     loce =  wrf_user_latlon_to_ij(f, ip_latse, ip_lonse)
 
     x_start = locs(1)
     x_end   = loce(1)
     y_start = locs(0)
     y_end   = loce(0)
 
     map=wrf_map_zoom(wks,f,mpres,y_start,y_end,x_start,x_end)

;  add province
   resp  = ARWres
   resp@gsLineColor      = "black"                     ; color of lines
   resp@gsLineThicknessF = 1.0                       ; thickness of lines
   dum = new(210,"graphic")
;   dirdat = "/disk2/zcx/wrf2/CHINA_mp/NEW/"
;   colorarr = random_uniform(1,150,204)
   do i = 1, 210
    yc   = asciiread(dirdat + i + ".dat",-1,"float")
    np = dimsizes(yc)
    m  = np/2
    xlat = new(m,float)
    xlon = new(m,float)
    xlat(0:m-1) = yc(0:np-2:2)
    xlon(0:m-1) = yc(1:np-1:2)
;   print(np)
;   resp@gsFillColor = colorarr(i-1)
    if(i.ge.204) then  ; adding Yellow River and Yangtz River
     resp@gsLineColor      = "chocolate"
     resp@gsLineThicknessF = 2.0
     resp@gsLineDashPattern = 0
    end if
    dum(i-1) = gsn_add_polyline(wks,map,xlon,xlat,resp)
;   dum(i-1) = gsn_add_polygon(wks,map,xlon,xlat,resp)
    delete(xlat)
    delete(xlon)
    delete(yc)
    delete(np)
   end do
;   draw(map)
;   frame(wks)
;   end province
; create colors
;*************************************************
  colors = (/ (/255,255,255/),(/0,0,0/),(/255,255,255/), (/166,242,143/), (/61,186,61/), (/97,184,255/), (/0,0,225/), (/250,0,250/), (/128,0,64/) /) * 1.0 ; we multiply by 1 to make colors float

  colors = colors/255.               ; normalize (required by NCL)

  gsn_define_colormap(wks, colors)   ; generate new color map
;***************************************************
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        
; What times and how many time steps are in the data set?
        
  ntimes = nfs      ; number of times in the file
  Times = f->Times
  Time_s = chartostring( Times )
   print("Time_s=" + Time_s)
  time_xa= wrf_times_c(Times,3)
   print("time_xa=" + time_xa)
  ac(ifs) = gmt2local(time_xa)  ; convert into Beijing Times
    ARWres@TimeLabel = ac(ifs)   ; Set Valid time to use on plots
;    ARWres@TimeLabel = Time_s   ; Set Valid time to use on plots
    print("ac(ifs)=" + ac(ifs))
    print("ARWres@TimeLabel = " + ARWres@TimeLabel )
    
  Times2 = f2->Times
  Time_s2 = chartostring( Times2 )
   print("Time_s2=" + Time_s2)
  time_xa2= wrf_times_c(Times2,3)
   print("time_xa2=" + time_xa2)
  ac2(ifs) = gmt2local(time_xa2)  ; convert into Beijing Times
    print("ac2(ifs)=" + ac2(ifs))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
  ; Get sea surface pressure, surface wind and surface temperature
    rain_exp = wrf_user_getvar(f,"RAINNC",0)-wrf_user_getvar(f2,"RAINNC",0)
    rain_con = wrf_user_getvar(f,"RAINC",0)-wrf_user_getvar(f2,"RAINC",0)
    rain_tot_24h = rain_exp + rain_con
    rain_tot_24h@description = "24-hour Total Precipitation : " + ac2(ifs) + " TO " + ac(ifs)
         ; Plotting options for Precipitation
        opts_r = ARWres
        opts_r@UnitLabel            = "mm"
        opts_r@cnLevelSelectionMode = "ExplicitLevels"
        opts_r@cnLevels             = (/ 0.1, 10., 25., 50., 100., 250./)
        opts_r@cnInfoLabelOn        = False
        opts_r@cnConstFLabelOn      = False
        opts_r@cnFillOn             = True
        opts_r@gsnFrame             = False
        opts_r@gsnMaximize          = True
        opts_r@gsnPaperOrientation = "portrait" 
        contour_tot_24h = wrf_contour(a,wks, rain_tot_24h(y_start:y_end,x_start:x_end), opts_r)
        delete(opts_r)
 ; MAKE PLOTS
        wrf_map_overlay(wks,map,contour_tot_24h,True)
     
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
  end do        ; END OF TIME LOOP
end