Dear ncl users,
I calculate the windpop variable for each grid point using do loops..
In order to plot the windpop values, I need to convert the windpop from 1-d
to 3-d (nlev,nlat,mlon)..
I use the onedtond function and the new variable wp has the right dimensions
(1x89x184). However, the values did not correspond with the right values..
The first value which belongs to first grid point (0,0,0) was repeated for
16376 times (up to (0,89,184) grid point), then the same is applied to the
second value etc.
The desired results would be (0,0,0) 1st value, (0,1,0) 2nd etc
I can't understand where is the error..
Has anybody got any idea?
Any help would be appreciated
The source files are available at
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1961-197
0_wss.nc.gz
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1971-198
0_wss.nc.gz
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1981-199
0_wss.nc.gz
for free download.
Attached script
;*************************************************
; Plot wind power potential over Mediterranean
;*************************************************
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
;=======================================================
; Variables definition and files handling
;=======================================================
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
pollat = 39.25
pollon = -162.0
;======================================================
; Choose time period
;======================================================
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
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
dimx = dimsizes(ws1)
nlev = dimx(1)
nlat = dimx(2)
mlon = dimx(3)
; npts = nlat*mlon
; wpower = new((/nlev,nlat,mlon/),"string")
npt = -1
do kl = 0,nlev-1
do nl = 0,nlat-1
do ml = 0,mlon-1
data = ws1(:,kl,nl,ml)
;======================================================================
; Calculate probability function
;======================================================================
opt = True
opt@bin_spacing = 0.2
max_ws1 = max(data)
min_ws1 = 0.
nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number
wspdf = pdfx(data,nbins,opt)
; printVarSummary(wspdf)
;=======================================================================
; Calculate cumulative density function
;=======================================================================
xcumsum = cumsum(wspdf,2)
xcumsum_Div = 99.99999 ; in order to avoid
xcumsum = where(xcumsum.ne.100,xcumsum,xcumsum_Div) ; impossible division
;=======================================================================
; Calculate Weibull parameters using least square method
;=======================================================================
yi = log(log(1.0/(1.0 - (xcumsum/100))))
if(any(isnan_ieee(yi))) then ; in order to avoid nan result
value = 2.230200812934578 ; in yi estimation
replace_ieeenan(yi,value,0)
yi@_FillValue = value
end if
yi_Zero = 0
yi = where(yi.ne."-inf",yi,yi_Zero) ; set missing (zero) value -inf
umaxiall = (/0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,\
3.,3.2,3.4,3.6,3.8,4.,4.2,4.4,4.6,4.8,5.,5.2,5.4,5.6,5.8,\
6.,6.2,6.4,6.6,6.8,7.,7.2,7.4,7.6,7.8,8.,8.2,8.4,8.6,8.8,\
9.,9.2,9.4,9.6,9.8,10.,10.2,10.4,10.6,10.8,11.,11.2,11.4,\
11.6,11.8,12.,12.2,12.4,12.6,12.8,13.,13.2,13.4,13.6,13.8,\
14.,14.2,14.4,14.6,14.8,15.,15.2,15.4,15.6,15.8,16.,16.2,\
16.4,16.6,16.8,17.,17.2,17.4,17.6,17.8,18.,18.2,18.4,18.6,\
18.8,19.,19.2,19.4,19.6,19.8,20.,20.2,20.4,20.6,20.8,21.,\
21.2,21.4,21.6,21.8,22.,22.2,22.4,22.6,22.8,23.,23.2,23.4,\
23.6,23.8,24.,24.2,24.4,24.6,24.8,25.,25.2,25.4,25.6,25.8,\
26.,26.2,26.4,26.6,26.8,27.,27.2,27.4,27.6,27.8,28.,28.2,28.4,\
28.6,28.8,29.,29.2,29.4,29.6,29.8,30.,30.2,30.4,30.6,30.8,\
31.,31.2,31.4,31.6,31.8,32.,32.2,32.4,32.6,32.8,33.,33.2,\
33.4,33.6,33.8,34.,34.2,34.4,34.6,34.8,35.,35.2,35.4,35.6,\
35.8,36.,36.2,36.4,36.6,36.8,37.,37.2,37.4,37.6,37.8,38.,\
38.2,38.4,38.6,38.8,39.,39.2,39.4,39.6,39.8,40./)
umaxi = umaxiall(0:nbins-1)
xi = log(umaxi)
xiyi= xi*yi
sumxy = sum(xiyi(0:nbins-1))
sumy = sum(yi(0:nbins-1))
sumx = sum(xi(0:nbins-1))
xi2 = xi*xi
sumx2 = sum(xi2(0:nbins-1))
b1 = sumy * sumx2
b2 = sumx * sumxy
b3 = nbins * sumx2
b4 = sumx * sumx
b = (b1 - b2) / (b3 - b4)
k1 = nbins * sumxy
k2 = sumx * sumy
k3 = nbins * sumx2
k4 = b4
;========================================================================
; Weibull Shape Parameter k
;========================================================================
k = (k1 - k2) / (k3 - k4)
;========================================================================
; Weibull Scale Parameter c
;========================================================================
ratio = - (b / k)
c = exp(ratio)
;========================================================================
; Calculate wind speed probability density function from 5 up to 25 m/s
;========================================================================
df1 = (k / c)
if(nbins.le.125) then
df2 = (umaxiall(25:nbins) / c) ^ (k-1)
df3 = -((umaxiall(25:nbins) / c) ^ k)
else
df2 = (umaxiall(25:125) / c) ^ (k-1)
df3 = -((umaxiall(25:125) / c) ^ k)
end if
df4 = exp(df3)
df = df1 * df2 * df4
numdf = (dimsizes(df)) - 1
;========================================================================
; Calculate Wind Power Potential
;========================================================================
ginol5a = umaxiall(25) ^ 3 ; epeidh epilegw thn akrh toy bin
if(nbins.le.125) then
ginol25a = umaxiall(nbins) ^ 3 ; to oloklirwma ypologizetai basei
else
ginol25a = umaxiall(125) ^ 3
end if
ginol5b = df(0) * 0.2 ; methodoy toy Simpson
ginol25b = df(numdf) * 0.2
ginol5 = ginol5a * ginol5b
ginol25 = ginol25a * ginol25b
if(nbins.le.125) then
ginol1 = (umaxiall(25:nbins))^3
else
ginol1 = (umaxiall(25:125))^3
end if
ginol2 = df * 0.2
ginol = ginol1 * ginol2
sumol0 = sum(ginol)
sumol = sumol0 + (0.5 * ginol25) - ginol5
windpop = 0.5 * 1.225 * sumol
; npt = npt+1
wp = onedtond(windpop,(/nlev,nlat,mlon/))
; printVarSummary(wp)
delete([/wspdf,data,xcumsum,umaxi,yi,xi,xiyi,xi2,df2,df3,df4,\
df,ginol1,ginol2,ginol/])
end do
end do
end do
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 : Thu Jul 25 2013 - 21:02:42 MDT