[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