Dear ncl users,
I would like to calculate the pdf of wind speed for a grid node.
When I calculated it for a grid point everything went well.
The number of bins, and the maximum wind speed for each grid point have been
estimated correctly (line number 92, 89 respectively), for each grid point
throughout the time period (10957 records)..
Moreover, the printVarSummary for max_ws1 and nbins give me the dimensions..
Number of dimensions 3 (1 x 19 x 17) -> (height x rlat x rlon)
As you can see at the attached script, I tried to solve the problem using do
loops but unfortunately without success.
Has anyone got any idea, how to calculate the pdf for each grid point..
Any help would be appreciated..
Ioannis
;*************************************************
; PDF Calculation
;*************************************************
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"
begin
DataDir = "/data3/koletsis/Current/A2/"
Files = systemfunc ( " ls " + DataDir + \
"C4IRCA3_A2_ECHAM5_DM_*.nc")
numFiles = dimsizes(Files)
print(Files)
print(" ")
f1 = addfiles(Files,"r")
time = f1[:]->time
ws = f1[:]->wss
ws@_FillValue = 1.e+30
lat2d = f1[0]->lat
lon2d = f1[0]->lon
time2 = ut_calendar(time, 0)
time3 = ut_calendar(time, -2) ; Date format
year = tointeger(time2(:,0)) ; Convert to integer
month = tointeger(time2(:,1))
day = tointeger(time2(:,2))
beg_year = 1961
end_year = 1990
beg_month = 01
end_month = 12
beg_day = 01
end_day = 31
ymdStrt = beg_year*10000+beg_month*100+beg_day
ymdLast = end_year*10000+end_month*100+end_day
print("ymdStrt="+ymdStrt+" ymdLast="+ymdLast)
itime = ind(time3.ge.ymdStrt.and.time3.le.ymdLast)
if (ismissing(itime(0))) then ; Error message
print("itime is missing:something is going wrong")
exit
end if
;======================================================================
; Select the subregion (Mediterranean Area)
;======================================================================
; latS = 30
; latN = 46
; lonW = -10
; lonE = 41
latS = 36
latN = 40
lonW = 22
lonE = 26
ji = region_ind(lat2d,lon2d,latS,latN,lonW,lonE)
jStrt = ji(0)
jLast = ji(1)
iStrt = ji(2)
iLast = ji(3)
llat2d = lat2d(jStrt:jLast,iStrt:iLast)
llon2d = lon2d(jStrt:jLast,iStrt:iLast)
ws2 = ws(itime,:,jStrt:jLast,iStrt:iLast)
ws1 = 1.2697*ws2 ; wind speed at hub height
; print(ws1)
; printVarSummary(ws2)
; printVarSummary(ws1)
;======================================================================
; Calculate probability function
;======================================================================
opt = True
opt@bin_spacing = 0.2
max_ws1 = dim_max_n(ws1,0)
min_ws1 = 0.
nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number
; print(nbins)
dimx = dimsizes(ws1)
nlev = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
do kk = 1, 10957
do nl = 0, nlat-1
do ml = 0, mlon-1
do nt = 0, nlev-1
wspdf=(pdfx(ws1(kk,nt,nl,ml),nbins,opt))
end do
end do
end do
end do
print(wspdf)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Fri Jul 12 2013 - 16:37:39 MDT