[Wrf-users] Fwd: Data analysis using NCL

Sara Fenech sarafenech at gmail.com
Wed Jul 2 17:21:10 MDT 2014


---------- Forwarded message ----------
From: Sara Fenech <sarafenech at gmail.com>
Date: Wed, Jul 2, 2014 at 12:44 PM
Subject: Re: [Wrf-users] Data analysis using NCL
To: Li Xianxiang <lixx at smart.mit.edu>
Cc: "wrf-users at ucar.edu" <wrf-users at ucar.edu>


Hi Li

First of all thanks for your reply. I have altered my script (shown below)
and it seems to work fine when opening the wrfout file for domain 1 (48km).
When I tried running the script for wrfout file of domain 2 (9,6km nest) I
got very unrealistic huge values (billions m/s).

Do you have any idea what the problem could be?

Kind Regards
Sara

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

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/wrf/WRFUserARW.ncl"

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


begin
;********************************************************
; read in netCDF file and make a loop for all time steps
;********************************************************
  in     = addfile("/home/wrf/OUTPUT/run_PBL2winter/wrfout_d01_2010-12-01.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

;T100 = new(ntimes,float)            ; creation of a T vector at each time
step
;P98 = new(ntimes,float)               ; creation of a P vector at each
time step
;Q = new(ntimes,float)               ; creation of a Q vector at each time
step
;rh100 = new(ntimes,float)               ; creation of a RH vector at each
time step
windspd100 = new(ntimes,float) ; creation of a Windspeed vector at each
time step
;print(ntimes)

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 = 35.99725
Longitude = 14.36775

res = True
res at 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)

 ;print("X location is: " + x)            ; print the value of X at the
screen
 ;print("Y location is: " + y)            ; print the value of Y at the
screen



;*************************************************************************************
;  - extract wind, Temperature, Pressure, relative humidity and height
coordinates-  *
;*************************************************************************************

; 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

; Temperature

 ;T_in = wrf_user_getvar(in, "tc", time)      ; Extract sea level pressure
(°C)

; Surface Pressure

 ;P_in = wrf_user_getvar(in, "pressure", time)      ; Extract sea level
pressure (Pa)



; Relative humidity

 ;rh_in = wrf_user_getvar(in, "rh", time)      ; Extract relative humidity
(%)

;print(rh_in)


;*******************************************************************************
;     - Interpolate wind speed and wind direction at 100m height -
 *
;*******************************************************************************



     ; Interpolate U,V to 100 Meters
      u_plane  = wrf_user_intrp3d( u,height,"h", 100,0.,False)
      v_plane  = wrf_user_intrp3d( v,height,"h", 100,0.,False)

      ; Calculate Wind Speed from Vectors
      spd = (u_plane*u_plane + v_plane*v_plane)^(0.5)
      windspd100(it)=spd(x,y)

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

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

;************************************************************
;  - Interpolate temperature at 100m height -                *
;************************************************************
  ; T_plane  = wrf_user_intrp3d( T_in,height,"h", 98,0.,False)
  ; T100(it)=T_plane(x,y)
;************************************************************
;  - Interpolate P at 98m height -                       *
;************************************************************
  ; P_plane  = wrf_user_intrp3d( P_in,height,"h", 98,0.,False)
  ; P98(it)=P_plane(x,y)
;************************************************************
;  - Interpolate rh at 98m height -                       *
;************************************************************

  ; rh_plane  = wrf_user_intrp3d( rh_in,height,"h", 98,0.,False)
  ; rh100(it)=rh_plane(x,y)





end do



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


print("  Time       Wind_speed_100m        ")

do it = 0,ntimes-1


 print (sprintf("%5.0f",it)    +" " \
                   +sprintf("%23.2f", windspd100(it)) +"  " )





end do                     ; end of time loop

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


    ;print(windspd100(it))

    line = new(721,string)   ; new(ntimes+1,string)    ; declare a string
array

    ; print headers
    line(0) = "  Time     Wind_speed        "

    do it = 0,ntimes-1

    line(it+1) =    sprintf("%5.0f",it)    +" " \
                  +sprintf("%23.2f", windspd100(it)) +"  "

    end do                     ; end of time loop

    print(""+line)         ; print output to console
    asciiwrite("PBL2AhraxDec.txt",line)      ; write output to file



end




On Mon, Jun 30, 2014 at 5:59 PM, Li Xianxiang <lixx at smart.mit.edu> wrote:

> Hi Sara,
>
> Basically your error is to interpolate with a single value (rather than an
> array). Here
>
> ua    = ua_in(0:0,y,x)
>
> is just a scalar rather than an array. It should be
>
> ua = us_in(:, y,x)
>
> The same applies to other variables. You can use printVarSummary to check
> the dimension of a variable.
>
> Other thoughts:
> 1. wrf_ll_to_ij returns the index of the point in the domain, starting
> from 1. But in ncl index of an array starts from 0. So usually we should
> subtract 1 from the i and j. However, I may be wrong here. You'd better
> check.
>
> 2. I remember wrf_user_getvar(nc_file, "ua", time) returns the unstaggered
> u. You can check.
>
> Xianxiang
>
>
> On 30 Jun, 2014, at 11:38 PM, "Sara Fenech" <sarafenech at gmail.com> wrote:
>
> Hi all
>
> I am very new to NCL and I've been trying to write a simple script but I'm
> having some problems. Well I have done 7 runs using wrf having hourly data
> for 3 months. Now I would like to extract the wind speed for each time step:
> 1) at a particular Lat and Lon (I managed this)
> 2) at a predefined height (in m)
>
> Now I have managed to plot the wind speed for all my domain at a
> particular and I have managed to get wind speeds for a particular lat and
> lon but I cannot manage to obtain both at once.
>
> Thus to sum it up I would like one value for wind speed for each time step
> (where lat lon and lev are predefined.)
>
> The last thing I would like to ask is how I can save the extracted data
> into a file.
>
> Thanks a lot for all you help. I'm attaching the script which is not
> working maybe you can spot what's wrong (Something to do with dimensions)
> together with the error message at the bottom.
>
> Thanks Again
>
> Sara
>
>
> ;******************************************************************************
> ;*NCL script to generate netCDF file containing wind speed at specific
> location    *
> ;* Example for wrfout_d01_2010-12-01.nc
>  *
> ;*                                        *
>
> ;************************************************************************************
>
>
> ; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------
>
> 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/wrf/WRFUserARW.ncl"
>
> ; --------------  BEGINING OF NCL SCRIPT ----------------
>
>
> begin
> ;********************************************************
> ; read in netCDF file and make a loop for all time steps
> ;********************************************************
>
>   a = addfile("/home/wrf/OUTPUT/run_PBL2winter/wrfout_d01_2010-12-01.nc
> ","r")
>
>
> times = wrf_user_list_times(a)            ; get times in the file
> ntimes = dimsizes(times)            ; number of times in the file
> wind_speed = new(ntimes,float) ; creation of a windspeed vector at each
> time step
> ;print(ntimes)
>
> do it = 0,ntimes-1                  ;Loop for the time: it= starting time
> print("Working on time " + it )
> time = it
>
>
> ;************************************************
> ;  - Select lon & lat of the point of interest -
> ;************************************************
>
>
>  res = True
>  res at returnInt = True  ; False : return real values, True: return
> interger values
>  lat = 42.47                  ; Latitude of the point of interest
>  lon = 12.98                  ; Longitude of the point of interest
>  point = wrf_user_ll_to_ij(a,lon,lat,res)       ;
> wrf_user_ll_to_ij(nc_file,lon,lat,opt)
>
>  x = point(1)
>  y = point(0)
>
>  print("X location is: " + x)            ; print the value of X at the
> screen
>  print("Y location is: " + y)            ; print the value of Y at the
> screen
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; The specific height levels that we want the data interpolated to.
>   height_levels = (/ 1039./)   ; height levels to plot - in meter
>   nlevels       = 1     ; number of height levels
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
>
> ;************************************************
> ;  - extract wind components and unstagger them -
> ;************************************************
>
>
> ; ----- Hence U and V are staggered, we need to  average them to get the
> point on the mass grid using wrf_user_unstagger
>
>   U     = wrf_user_getvar(a, "U", time)
>   ua_in = wrf_user_unstagger(U,U at stagger)
>   ua    = ua_in(0:0,y,x)             ;ua(bottom_up,grid_lon,grid_lat)
>   V     = wrf_user_getvar(a, "V", time)
>   va_in = wrf_user_unstagger(V,V at stagger)
>   va    = va_in(0:0,y,x)            ;va(bottom_up,grid_lon,grid_lat)
>   W     = wrf_user_getvar(a, "W", time)
>   z_in = wrf_user_unstagger(W,W at stagger)
>   z    = z_in(0:0,x,y)             ;z(bottom_up,grid_lon,grid_lat)
>
>   wind_speed(it) = sqrt(ua^2+va^2)
>   ;copy_VarCoords(ua,wind_speed(it))     ; copy coord vars to speed
>   ;print(wind_speed(it))
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;Interpolating to the requested height
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>       u_plane  = wrf_user_intrp4d( ua,z,"h",height_levels,0.,False)
>       v_plane  = wrf_user_intrp3d( va,z,"h",height_Levels,0.,False)
>       spd_plane =wrf_user_intrp3d
> (wind_speed(it),z,"h",height_Levels,0.,False) ; speed in m/sec
>
>    print(spd_plane(it))
>
> end do                     ; end of time loop
> ;************************************************
> ;  - Write wind speed in ascii file -
> ;************************************************
>
>
>
>
>  ; XLAT = lat              ; (nlat)
>  ; YLON = lon              ; (mlon)
>   ;windspeed   = wind_speed             ; (ntim,nlat,mlon)
>
>   ;dimwindspeed  = dimsizes(windspeed)
>   ;ntim  = dimwindspeed(0)
>   ;nXLAT  = dimwindspeed(1)
>   ;mYLON  = dimwindspeed(2)
>
>   ;npts  = nXLAT*mYLON        ; total number of grid points
>
>   ;fName = "foo.txt"
>   ;data  = new( npts, "string")
>
>   ;npt   = -1
>
>   ;do nl=0,nXLAT-1
>    ; do ml=0,mYLON-1
>
>     ;   npt  = npt + 1
>      ;  data(npt) = sprinti("%0.5i", (npt+1) )
>       ; data(npt) = data(npt) + sprintf("%7.1f ",XLAT(nl))
>        ;data(npt) = data(npt) + sprintf("%7.1f ",YLON(ml))
>
>       ;do nt=0,ntim-1
>        ;  data(npt) = data(npt) + sprintf("%10.3f ", windspeed(nt,nl,ml))
>       ;end do
>     ;end do
>   ;end do
>
>   ;asciiwrite (fName , data)
>
>
> ; --------------  END OF NCL SCRIPT ----------------
>
> end
>
> *Error Message*
>
> fatal:Subscript out of range, error in subscript #0
> fatal:An error occurred reading dims
> fatal:["Execute.c":8567]:Execute: Error occurred at or near line 223 in
> file $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl
>
> fatal:["Execute.c":8567]:Execute: Error occurred at or near line 84 in
> file SpLev.ncl
>
>
> _______________________________________________
> Wrf-users mailing list
> Wrf-users at ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/wrf-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/wrf-users/attachments/20140703/9ac038ce/attachment-0001.html 


More information about the Wrf-users mailing list