;************************************************
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/csm/shea_util.ncl"  
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;************************************************

  a0 = addfile("/home/swati/work/analysis/2_cities/Patna/14aug2014/wrf_out_data/3rd_domain/13_aug.nc","r")
  a1 = addfile("/home/swati/work/analysis/2_cities/Patna/14aug2014/wrf_out_data/3rd_domain/14_aug.nc","r")
  a2 = addfile("/home/swati/work/analysis/2_cities/Patna/14aug2014/wrf_out_data/3rd_domain/15_aug.nc","r")

; The specific pressure levels that we want the data interpolated to.
  pressure = (/ 850 /)   ; pressure levels to plot
    p  = wrf_user_getvar(a0,"pressure",-1) ; used to mention vertical coordinate to interplate to

; First get the variables we will need        
;----------------------------------    a0   ----------------------------------------
    u0  = wrf_user_getvar(a0,"ua",-1)        ; u averaged to mass points
    v0  = wrf_user_getvar(a0,"va",-1)        ; v averaged to mass points

    uu0  = wrf_user_intrp3d( u0,p,"h",pressure,0.,False)
    vv0  = wrf_user_intrp3d( v0,p,"h",pressure,0.,False)

    u_plane0  = dim_avg_Wrap(uu0(south_north|:,west_east|:, Time|:) )
    v_plane0  = dim_avg_Wrap(vv0(south_north|:,west_east|:, Time|:) )
    spd0     = (u_plane0*u_plane0 + v_plane0*v_plane0)^(0.5)

;----------------------------------    a1   ----------------------------------------
    u1  = wrf_user_getvar(a1,"ua",-1)        ; u averaged to mass points
    v1  = wrf_user_getvar(a1,"va",-1)        ; v averaged to mass points

    uu1  = wrf_user_intrp3d( u1,p,"h",pressure,0.,False)
    vv1  = wrf_user_intrp3d( v1,p,"h",pressure,0.,False)

    u_plane1  = dim_avg_Wrap(uu1(south_north|:,west_east|:, Time|:) )
    v_plane1  = dim_avg_Wrap(vv1(south_north|:,west_east|:, Time|:) )
    spd1     = (u_plane1*u_plane1 + v_plane1*v_plane1)^(0.5)

;----------------------------------    a2   ----------------------------------------

    u2  = wrf_user_getvar(a2,"ua",-1)        ; u averaged to mass points
    v2  = wrf_user_getvar(a2,"va",-1)        ; v averaged to mass points

    uu2  = wrf_user_intrp3d( u2,p,"h",pressure,0.,False)
    vv2  = wrf_user_intrp3d( v2,p,"h",pressure,0.,False)

    u_plane2  = dim_avg_Wrap(uu2(south_north|:,west_east|:, Time|:) )
    v_plane2  = dim_avg_Wrap(vv2(south_north|:,west_east|:, Time|:) )
    spd2     = (u_plane2*u_plane2 + v_plane2*v_plane2)^(0.5)