load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

begin

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;WRF PORTION;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 
a= addfile("/media/mughal/omer data disk/WRf_RUNS/ECMWF_full_month/ECMWF_final_1/wrfout_d04_2009-07-10_00:00:00.nc","r") 
  
wks = gsn_open_wks("x11","gsun08n") 

Times = a->Times                     ;reading in time variable from the wrfout file
Time_3 = wrf_times_c( Times, 3 )     ; opt=3 returns units of "yyyymmddhh" (type "integer")
                                     ;Read the character "Times" variable from a netCDF file. Convert to 					    
                                     ;different units as specified by opt.  
time=Time_3(144)                     ;constraining the value to be read into just one day at the first hour
print("Working on time: " + time )

 
time2=144                           ;values to be read at the specified time  
;;;;Reading in the required variables;;;;;;;;;;

u= wrf_user_getvar(a,"ua",time2)  
v= wrf_user_getvar(a,"va",time2)  
spd = (u*u + v*v)^(0.5) 
;u = u*1.94386                       ; winds now in kts
      ;v = v*1.94386                       ; winds now in kts
      ;u@units = "kts"
      ;v@units = "kts"
;;;;Assigining resources;;;;;;;;;;;;;
res            =True
res@TimeLabel = time 
res@MainTitle = "REAL-TIME WRF"
res@TimeLabel = time 
;res@gsnDraw   = False
;res@gsnFrame  = False
pltres = True
pltres@NoTitles =True
mpres  = True
mpres@mpGeophysicalLineColor = "Black"
mpres@mpNationalLineColor    = "Black"
mpres@mpUSStateLineColor     = "Black"



;res@ContourParameters = (/ 15., 60., 5. /)
res@cnFillOn          = True
contour = wrf_contour(a,wks,spd(2,:,:),res) 
    
opts=True
opts@FieldTitle        = "Wind"       ; Overwrite Field Title
opts@NumVectors        = 47           ; wind barb density

vector =  wrf_vector(a,wks,u(2,:,:),v(2,:,:),opts)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;LIDAR PORTION;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

b=addfile("./20090711_000239_01_g.nc","r")

lon=b->yp(0,0,:)
lat=b->xp(0,:,0)

u1=b->vrect_x(:,:,:)
 u1@_FillValue = -99.0


 v1=b->vrect_y(:,:,:)       
 v1@_FillValue = -99.0
 u1!0="level"
 u1!1="lon"
 u1!2="lat"  
 u1&lat=lat
 u1&lon=lon 
 u1&lat@units="degrees_north"
 u1&lon@units="degrees_east"

 v1!0="level"
 v1!1="lon"
 v1!2="lat"
 v1&lat=lat
 v1&lon=lat
 v1&lat@units="degrees_north" 
 v1&lon@units="degrees_east" 
;u = u*1.94386                       ; winds now in kts
      ;v = v*1.94386                       ; winds now in kts
      ;u@units = "kts"
      ;v@units = "kts"
vcres                 = True  
vcres@vfXCStartV  = lon              
vcres@vfXCEndV    = lon(dimsizes(lon)-1) 
vcres@vfYCStartV  = lat
vcres@vfYCEndV    = lat(dimsizes(lat)-1) 
vcres@tiXAxisString = "longitude"
vcres@tiYAxisString = "latitude"
 ;vcres@vcGlyphStyle = "WindBarb"
;vcres@vcMonoWindBarbColor     = False            ; color barbs by scalar 
 ; vcres@vcMinDistanceF          = 0.025            ; thin out windbarbs
vcres@vcMinFracLengthF = 0.33
vcres@vcRefMagnitudeF  =80.0
vcres@vcRefLengthF     = 0.045
vcres@vcMonoLineArrowColor  = False
vcres@gsnDraw = False			
vcres@gsnFrame = False

vplot = gsn_vector(wks,u1(4,::10,::10),v1(4,::10,::10),vcres)
 

;plot = wrf_map_overlays(a,wks,(/contour, vector,vplot/),pltres,mpres)

plot = wrf_map_overlays(a,wks,(/contour, vector/),pltres,mpres)

overlay(plot,vplot)


draw(plot)

frame(wks)  

maximize_output(wks,True)
end