;************************************************* ; Time-Averaged and Plot wind speed fields ;************************************************* 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 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) ; ws2 = ws(itime,:,40,145) ws1 = 1.2697*ws2 ; wind speed at hub height ; print(ws1) ; printVarSummary(ws1) ;====================================================================== ; Calculate probability function ;====================================================================== opt = True opt@bin_spacing = 0.2 ; opt@bin_min = 0. ; opt@bin_max = 25. ; max_ws1 = dim_max_n(ws1,0) ; max_ws1 = max(ws1,0) ; min_ws1 = 0. ; nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number ntim = dimsizes(ws(:,0,0,0)) nlev = dimsizes(ws(0,:,0,0)) nlat = dimsizes(llat2d(:,0)) mlon = dimsizes(llon2d(0,:)) do kl = 0,nlev-1 do nl = 0, nlat-1 do ml = 0, mlon-1 max_ws1 = dim_max_n(ws1,0) min_ws1 = 0. nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number pdf_grdpt = pdfx(ws1(:,kl,nl,ml),nbins,opt) printVarSummary(pdf_grdpt) end do end do end do ; wspdf=onedtond(pdfx(ws1,nbins,opt),(/nbins,nlat,nlon/)) ; print(max_ws1) ; printVarSummary(wspdf) end