begin
a = addfile("/scale_wlg_nobackup/filesets/nobackup/uoa02450/2017_Model Results/wrfout_d03_2017JFMAMJJASOND_00:00:00_Time.nc","r")
type = "pdf"

wks = gsn_open_wks(type,"Wind Power Density Jan-Dec 2017")

res = True
res@mpProjection  = "Mercator"    ; The default
res@mpDataBaseVersion = "HighRes"

times = wrf_user_getvar(a,"times",-1) 
ntimes = dimsizes(times)       
       
u = wrf_user_getvar(a,"ua",-1) 
v = wrf_user_getvar(a,"va",-1)   
spd = (u*u + v*v)^(0.5)

spd_avg = dim_avg_n(spd(0:365,:,:,:),0)

wpd = 0.5*1.225*((spd_avg)^(3))
wpd@description = "Wind Power Density"
wpd@units = "W m-2" 

res = True               
opts =  res
opts@cnFillOn = True  
opts@cnLinesOn = True
opts@ContourParameters = (/ 0., 500., 50/)
opts@gsnSpreadColorEnd = -3  
contour_wpd = wrf_contour(a,wks,wpd(0,:,:),opts)
delete(opts)

plot = wrf_map_overlays(a,wks,(/contour_wpd/),res,res)

end