load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"      ;Load the required library files
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.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/wrf/WRF_contributed.ncl" 

begin
begTime = get_cpu_time()
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	a = addfile("test.nc","r")
	xlat = a->latitude
	xlong = a->longitude
	nlat = dimsizes(xlat)
	nlon = dimsizes(xlong)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Extract single level parameters	
	tt2 = a->TMP_2maboveground
	t2 = tt2-273.15		;convert to C
	RH2 = a->RH_2maboveground
        
	t2@long_name = "Temperature @ 2m"
        t2@short_name = "T2"
        t2@units = "C"
        t2!0 = "time"
        t2!1 = "latitude"
        t2!2 = "longitude"
        RH2@long_name = "Relative Humidity @ 2m"
        RH2@short_name = "RH2"
        RH2@units = "%"
        RH2!0 = "time"
        RH2!1 = "latitude"
        RH2!2 = "longitude"

	lev = (/1000.,950.,925.,850.,700.,600.,500.,400.,300.,200.,100./)
	nlev = dimsizes(lev)
	t_1000 = a->TMP_1000mb
	t_950 = a->TMP_950mb
	t_925 = a->TMP_925mb
	t_850 = a->TMP_850mb
	t_700 = a->TMP_700mb
	t_600 = a->TMP_600mb
	t_500 = a->TMP_500mb
	t_400 = a->TMP_400mb
	t_300 = a->TMP_300mb
	t_200 = a->TMP_200mb
	t_100 = a->TMP_100mb
	ntime  = dimsizes(t_1000(:,0,0))
        t = new((/nlev,ntime,nlat,nlon/),"float")
	t = (/t_1000,t_950,t_925,t_850,t_700,t_600,t_500,t_400,t_300,t_200,t_100/)-273.15

	t@long_name = "Temperature"
        t@short_name = "tc"
        t@units = "C"
        t!0 = "level"
        t!1 = "time"
        t!2 = "latitude"
        t!3 = "longitude"


        u = new((/nlev,ntime,nlat,nlon/),"float")
	u_1000 = a->UGRD_1000mb
	u_950 = a->UGRD_950mb
	u_925 = a->UGRD_925mb
	u_850 = a->UGRD_850mb
	u_700 = a->UGRD_700mb
	u_600 = a->UGRD_600mb
	u_500 = a->UGRD_500mb
	u_400 = a->UGRD_400mb
	u_300 = a->UGRD_300mb
	u_200 = a->UGRD_200mb
	u_100 = a->UGRD_100mb
	u = (/u_1000,u_950,u_925,u_850,u_700,u_600,u_500,u_400,u_300,u_200,u_100/)
	
        u@long_name = "Zonal Wind"
        u@short_name = "ua"
        u@units = "m/s"
        u!0 = "level"
        u!1 = "time"
        u!2 = "latitude"
        u!3 = "longitude"

printVarSummary(u)
	v = new((/nlev,ntime,nlat,nlon/),"float")
	v_1000 = a->VGRD_1000mb
	v_950 = a->VGRD_950mb
	v_925 = a->VGRD_925mb
	v_850 = a->VGRD_850mb
	v_700 = a->VGRD_700mb
	v_600 = a->VGRD_600mb
	v_500 = a->VGRD_500mb
       	v_400 = a->VGRD_400mb
	v_300 = a->VGRD_300mb
	v_200 = a->VGRD_200mb
	v_100 = a->VGRD_100mb
	v = (/v_1000,v_950,v_925,v_850,v_700,v_600,v_500,v_400,v_300,v_200,v_100/)
	
        v@long_name = "Meridonial Wind"
        v@short_name = "va"
        v@units = "m/s"
        v!0 = "level"
        v!1 = "time"
        v!2 = "latitude"
        v!3 = "longitude"
printVarSummary(v)       


 	wd_1 = new((/nlev,ntime,nlat,nlon/),"float")  
 	ws_1 = new((/nlev,ntime,nlat,nlon/),"float")  	

	ws_1 = sqrt(u^2 +v^2)		
	wd_1 = wind_direction(u,v,1)	

	ws_1@long_name = "Wind Speed"
        ws_1@short_name = "ws"
        ws_1@units = "m/s"
        ws_1!0 = "level"
        ws_1!1 = "time"
        ws_1!2 = "latitude"
        ws_1!3 = "longitude"
	
	wd_1@long_name = "Wind Speed"
        wd_1@short_name = "wd"
        wd_1@units = "m/s"
        wd_1!0 = "level"
        wd_1!1 = "time"
        wd_1!2 = "latitude"
        wd_1!3 = "longitude"
        
        
        
        vvel = new((/nlev,ntime,nlat,nlon/),"float")
	vvel_1000 = a->VVEL_1000mb
	vvel_950 = a->VVEL_950mb
	vvel_925 = a->VVEL_925mb
	vvel_850 = a->VVEL_850mb
	vvel_700 = a->VVEL_700mb
	vvel_600 = a->VVEL_600mb
	vvel_500 = a->VVEL_500mb
	vvel_400 = a->VVEL_400mb
	vvel_300 = a->VVEL_300mb
	vvel_200 = a->VVEL_200mb
	vvel_100 = a->VVEL_100mb
	vvel = (/vvel_1000,vvel_950,vvel_925,vvel_850,vvel_700,vvel_600,vvel_500,vvel_400,vvel_300,vvel_200,vvel_100/)
	
	vvel@long_name = "Verical Wind"
        vvel@short_name = "vvel"
        vvel@units = "Pa/s"
        vvel!0 = "level"
        vvel!1 = "time"
        vvel!2 = "latitude"
        vvel!3 = "longitude"
printVarSummary(vvel)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Start the conversion of Pa/s to m/s for vertical wind
;omega_to_w(vvel,p,t)  "NCL Function" where p, t,vvel are in same dimension(t is in K, p is in Pa and vvel is in Pa/s)
       t       = t + 273.15              ; degC -> degK
       t@units = "degK"
                                         ; make p[*] conform to vvel[*][*][*][*]
       P       = conform(vvel, lev, 0)   ; P(:,:,:,:)
       P       = P*100                   ; hPa -> Pa
       P@units = "Pa"
       
       w       = omega_to_w(vvel, P, t)      ; w[*][*][*][*] (m/s)

       delete(P)                          ; no longer needed
;End the conversion
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

        w@long_name = "Verical Wind"
        w@short_name = "wa"
        w@units = "m/s"
        w!0 = "level"
        w!1 = "time"
        w!2 = "latitude"
        w!3 = "longitude"

printVarSummary(w)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
print("File creation and parameter writing starts")

	fw = addfile("gfsout_2015113000.nc", "c")
	tc = new((/nlev,ntime,nlat,nlon/),"float")
 
        ua = new((/nlev,ntime,nlat,nlon/),"float")
        va = new((/nlev,ntime,nlat,nlon/),"float")
        wa = new((/nlev,ntime,nlat,nlon/),"float")
        ws = new((/nlev,ntime,nlat,nlon/),"float")
        wd = new((/nlev,ntime,nlat,nlon/),"float")
        
	T2 = new((/ntime,nlat,nlon/),"float")
        rh2 = new((/ntime,nlat,nlon/),"float")


	fw->ua = u
	fw->va = v
	fw->wa = w	
	fw->ws = ws_1
	fw->wd = wd_1
	fw->T2 = t2
	fw->rh2 = RH2
print(begTime)
end