undef("reshape") function reshape(x,dims) begin return(onedtond(ndtooned(x),dims)) end ;************************************************* ; Calculation of PDF ;************************************************* 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/" DataDir = "./" Files = systemfunc ( " ls " + DataDir + \ "C4IRCA3_A2_ECHAM5_DM_25km_1961-1970_wss.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 = 1970 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) ntim = dimsizes(ws(:,0,0,0)) nlat = dimsizes(llat2d(:,0)) nlon = dimsizes(llon2d(0,:)) ws3 = ws(itime,:,jStrt:jLast,iStrt:iLast) ws2 = reshape(ws3,(/ntim,nlat*nlon/)) ws1 = 1.2697*ws2 ; wind speed at hub height printVarSummary(ws1) printMinMax(ws1,0) ;====================================================================== ; 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 wspdf = new((/ntim,max(nbins)/),double,1e30) printVarSummary(nbins) ; (nlat*nlon) printVarSummary(ws1) ; ntim x (nlat*nlon) printVarSummary(wspdf) do i=0,nlat*nlon-1 tmp = pdfx(ws1(:,i),nbins(i),opt) if(i.eq.0) then copy_VarAtts(tmp,wspdf) wspdf@_FillValue = tmp@_FillValue end if wspdf(i,0:dimsizes(tmp)-1) = tmp delete(tmp) end do printVarSummary(wspdf) printMinMax(wspdf,0) end