; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------

    load "/usr/share/ncarg/nclscripts/csm/gsn_code.ncl"
    load "/usr/share/ncarg/nclscripts/csm/gsn_csm.ncl"
    load "/usr/share/ncarg/nclscripts/csm/contributed.ncl"
    load "/usr/share/ncarg/nclscripts/wrf/WRFUserARW.ncl"
   load "/usr/share/ncarg/nclscripts/contrib/cd_string.ncl"

    ; --------------  BEGINING OF NCL SCRIPT ----------------


    begin
    ;********************************************************
    ; read in netCDF file and make a loop for all time steps
    ;********************************************************
      in     = addfile("./wrfout_d01_2014-01-03_06:00:00.nc","r")


    ;********************************************************
    ; Process all the time steps
    ;********************************************************

    times = wrf_user_list_times(in)     	; get times in the file
    ntimes = dimsizes(times)            	; number of times in the file
   
    windspd0 = new(ntimes,float)            	; creation of a Windspeed vector at each time step
    dir0 = new(ntimes,float) 

   do it = 0,ntimes-1				;Loop for the time: it= starting time
   time = it
 


    ;************************************************
    ;  - Select lon & lat of the point of interest -
    ;************************************************
    ; - The function wrf_user_ll_to_ij find the nearest grid point for a specific lat and lon

  Latitude =   -8.551083
    Longitude = 125.659373
    res = True      
    res@returnInt = True        		               ; False : return real values, True: return interger values
    point = wrf_user_ll_to_ij(in,Longitude,Latitude,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)
     
     x = point(0)
     y = point(1)

    



    ;*************************************************************************************
    ;  - extract wind
    ;*************************************************************************************
    
      ; Wind and Height

        u  = wrf_user_getvar(in,"ua",time)	        ; u averaged to mass points
        v  = wrf_user_getvar(in,"va",time)      	; v averaged to mass points
        height  = wrf_user_getvar(in, "z",time) 	; height is our vertical coordinate
        ter = wrf_user_getvar(in, "ter",time) 		; model terrain height (HGT_M, HGT_U, HGT_V)
     
    ; Conform data to Terrain Height   
       nheight = conform(height,ter,(/1,2/)) 		; assuming height is a 3d array and ter is a 2d array
       height = height - nheight


  
 	plane = (/ 0., 10./)            ; height levels to plot - in meter
         
         ; Interpolate U,V to 0 Meters
         	 u_plane  = wrf_user_intrp3d( u,height,"v",plane,0.,False)
         	 v_plane  = wrf_user_intrp3d( v,height,"v",plane,0.,False)
     	
          ; Calculate Wind Speed from Vectors
        	  spd = (u_plane*u_plane + v_plane*v_plane)^(0.5)
         	 windspd0(it)=spd(x,y)

          ; Wind direction at 0 Meters
       		r2d = 45.0/atan(1.0)       			; conversion factor (radians to degrees)
       		dir = atan2(u_plane,v_plane) * r2d + 180
       		dir0(it) = dir(x,y)

          ; Wind Speed
          spd@description = "Wind Speed"
          spd@units = "m/s"
          u_plane@units = "m/s"


  
    end do


    ;************************************************************
    ;  - Print the variables at the screen -                    *
    ;************************************************************

    npts=ntimes
    fName = "ncl0-10m.txt"
    data  = new( npts, "string")
       
     
print("Time(hour)    Wind_Direction(º)           Wind_Speed(m/s)  ")
     
do it = 0,ntimes-1
 
        print (sprintf  ("%5.0f",it)         +" "      \
       	  +sprintf ("               %4.1f", dir0(it))      +"  " \   
          +sprintf("%23.2f", windspd0(it)) +"  " )


    end do                     ; end of time loop

 asciiwrite (fName , data)

end