;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "WRFOptions.ncl"     ; set basic plot options here
;load "gsn_code.ncl"
;load "/home/cust100021_vol1/barry//barry/WRF_NCL/WRFPlot.ncl"
;load "WRFUserARW.ncl"
;load "SkewTFunc.ncl"
 begin

;   mpdatabaseversion= "Earth..4"

a = addfile("___FILE___.nc","r")

print(a)
; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"
;

; find ntimes
  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file
  print ("ntimes = " + ntimes)
;  ntimes=264
; ntimes=8
; create some variables
; ni-1,nj-1
  n_hours=24
  n_inc = 1
  inc = 1

  ter = wrf_user_getvar(a,"HGT",0)
  ter = ter - 10
 printMinMax(ter,0)
  dims=dimsizes(ter)
  print("dims(0) " + dims(0))
  print("dims(1) " + dims(1))
; u_arr=new((/dims(1),dims(0),n_hours/),float)
; v_arr=new((/dims(1),dims(0),n_hours/),float)
; a_arr=new((/dims(1),dims(0)/),float)
; b_arr;new((/dims(1),dims(0)/),float)
  coef = 0.25

;
; wtype = "png"
  wtype = "pdf"
; for png --> gives "error" for pdf plot
  wtype@wkWidth  = 2500    ; produces a better looking PNG
  wtype@wkHeight = 2500
  wks = gsn_open_wks(wtype,"wrf_israel_surf_12km.pdf")
; set up city labels
; cities = (/" Tel-Aviv"," Jerusalem"," Haifa"," Hermon","Arad","Eilat"," Med. Sea","Amman","Damascus","Nahariya","Dafna"/)
  cities = (/"  T-A","  J","  H","  H","  A","  E"," M-S","  A","  D-cus","  S",\
              "  E","  M-R", "  H", " A","  E-G","  N","  B-7","  E-A","  A", \
              "  N","  T","  Z","  A"/)
  lat    = (/    32.0670, 31.78300,32.808,33.368,31.2475,29.55,33.75,31.9454,33.4,32.967, \
                 31.65,30.617,31.53,31.677,31.450,32.334,31.25, 31.13,32.92,33.006,33.27,32.57, \
                 31.80/)
  lon    = (/    34.767,35.223,35.057,35.783,35.203515,34.95,33.75,35.9284,36.5,35.493, \
                 35.133,34.80,35.10,34.583,35.383,34.858,34.80,33.80,35.10,35.095,35.20,34.943, \
                 34.655/)


  ncities = dimsizes(cities)

  textid  = new(ncities,graphic)   ; Create array to hold text objects
  markid  = new(ncities,graphic)   ; Create array to hold marker objects

  txres               = True       ; Set up text resources
  txres@txFontHeightF = 0.01
  txres@txJust        = "CenterLeft"

  mkres               = True        ; Set up marker resources
  mkres@gsMarkerIndex = 16          ; Filled dot
  mkres@gsMarkerSizeF = .003          ; Size of dot

; find min,max and "avetc2"
  tc_whole = wrf_user_getvar(a,"T2",-1)   ; T2 in Kelvin
  min_tc2 = min(tc_whole)
  max_tc2 = max(tc_whole)
  print("min_tc2 = " + min_tc2)
  print("max_tc2 = " + max_tc2)
  avetc2=floattointeger(avg(tc_whole) -273.16)
  rngtc2 = tointeger(0.5*(max_tc2 + min_tc2)+2 - 273.16)
  print("rngtc2 = " + rngtc2)

;
;  ntimes=1
   do it = 0,ntimes-1,4
;  do it = 0,2
;  do it = 0,1,1
    time=it
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

;
    u  = wrf_user_getvar(a,"ua",time)      ; u averaged to mass points
    v  = wrf_user_getvar(a,"va",time)      ; v averaged to mass points

    if (it .eq. 0) then
      u10 = u(0,:,:)                       ; Use lowest level at time 0
      v10 = v(0,:,:)
    else
      u10 = wrf_user_getvar(a,"U10",time)  ; u at 10 m, mass point
      v10 = wrf_user_getvar(a,"V10",time)  ; v at 10 m, mass point
    end if
    u_plane = u10*3.6
    v_plane = v10*3.6
    tc  = wrf_user_getvar(a,"tc",time)     ; T in C
    tsk  = wrf_user_getvar(a,"TSK",time)     ; T in C
    td  = wrf_user_getvar(a,"td",time)     ; Td in C
    if ( it .eq. 0 ) then
      tc2 = tc(0,:,:)                      ; Use lowest T at time zero
    else
      tc2 = wrf_user_getvar(a,"T2",time)   ; T2 in Kelvin
      tc2 = tc2-273.16                     ; T2 in C
    end if
    tc_plane=tc2
;   tc_plane=tsk - 273.16
;   avetc2=floattointeger(avg(tc2))
    opts_tc = res
    opts_tc@FieldTitle         = "Surface Weather Map"
    opts_tc@UnitLabel          = "C"
   ; opts_tc@ContourParameters  = (/ -5., 5., 1./)
    opts_tc@ContourParameters  = (/ -rngtc2+avetc2, avetc2+rngtc2, 2./)
    opts_tc@cnFillOn           = True
    opts_tc@gsnSpreadColorEnd  = -3  ; End third from the last color in color map
;   opts_tc@NoHeaderFooter = True
    opts_tc@Footer = False
;   opts_tc@NoFooter = True
    contour_tc  = wrf_contour(a,wks,tc_plane,opts_tc)
;       Plotting options for SLP
       T  = wrf_user_getvar(a,"T",time) 
       th = T + 300.

       P  = wrf_user_getvar(a,"P",time) 
       PB  = wrf_user_getvar(a,"PB",time) 
       p  = ( P + PB )

       tk = wrf_tk( p , th )

       QVAPOR  = wrf_user_getvar(a,"QVAPOR",time) 
       PH  = wrf_user_getvar(a,"PH",time) 
       PHB  = wrf_user_getvar(a,"PHB",time) 
       var    = ( PH + PHB ) / 9.81
       dim    = dimsizes(var)
       z      = 0.5 * ( var(0:dim(0)-2,:,:) + var(1:dim(0)-1,:,:) )
       slvl   = wrf_slp( z, tk, p, QVAPOR )

    if(it.eq.0)then
    avepsl = floattointeger(avg(slvl))
    end if
    opts_psl = res
    res@gsnAddCyclic         = False   ; Don't add longitude cyclic pt.
    ;opts_psl@FieldTitle        = "Sea Level Pressure"
    ;opts_psl@UnitLabel         = "mb"
    opts_psl@ContourParameters = (/ -40+avepsl, avepsl+40, 4. /)
    opts_psl@cnLineColor       = "Grey"
    opts_psl@cnHighLabelsOn    = True
    opts_psl@cnLowLabelsOn     = True
   ;opts_psl@cnLineLabelBackgroundColor = -1
    opts_psl@cnLineLabelBackgroundColor = 0
    opts_psl@gsnContourLineThicknessesScale = 2.0
    opts_psl@NoHeaderFooter = True
    contour_psl = wrf_contour(a,wks,slvl,opts_psl)
      opts_psl@cnLineColor         = "NavyBlue"
      opts_psl@gsnContourLineThicknessesScale = 3.0


; Plotting options for Wind Vectors                 
      opts = res         
      opts@FieldTitle = "Wind"       ; overwrite Field Title
      opts@UnitLabel         = "km/hr"
      opts@NumVectors = 20           ; density of wind barbs
      opts@vcRefMagnitudeF = 5.0
      opts@vcRefLengthF=0.020
; turns off country boundaries
;     opts@mpOutlineBoundarySets = "NoBoundaries"
;   opts@mpOutlineBoundarySets = "Geophysical"
    opts@NoHeaderFooter = True
     vector      = wrf_vector(a,wks,u_plane, v_plane,opts)
;   plot = wrf_map_overlays(a,wks,(/contour_ter,vector/),False,False)
;First, in order to add markers and things to a WRF map plot, 
;you have to tell the wrf_map_overlays procedure to not remove the contour plot after it's done:
    pltres            = True   ; Set plot options
    pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove contours
    mpres             = True   ; Set map options (none set here)
  mpres@mpOutlineOn                 = True
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpUSStateLineThicknessF     = 3.0
  mpres@mpNationalLineThicknessF    = 3.0
  mpres@mpGeophysicalLineThicknessF = 3.0
  mpres@mpOutlineBoundarySets       = "AllBoundaries"


    plot = wrf_map_overlays(a,wks,(/contour_tc,vector,contour_psl/),pltres,mpres)
;   plot = wrf_map_overlays(a,wks,(/contour_tc,vector,contour_psl/),pltres,False)

;   markid = gsn_add_polymarker(wks,plot,lon,lat,mkres)
;   textid = gsn_add_text(wks,plot,cities,lon,lat,txres)


;Finally, in order to see the markers and the cities, you actually have to draw the map and advance the frame yourself.  The gsn_add_xxxx calls just attach the stuff to the map, but nothing gets drawn until you call:


    draw(plot)
    frame(wks)

      delete(opts)
   end do  ;       end of the time loop
end