; Concepts illustrated: ; ; - Calculating daily anomalies from the mean annual cycle ; - Getting the indices where data falls in a particular range ; - Converting "string" time values using cd_calendar ; - Writing data to a NetCDF file using the efficient method ; - Paneling six XY plots on the same page ; - Making all curves in an XY plot solid ; - Filling the areas of an XY curve above and below a reference line ; - Drawing a Y reference line in an XY plot ; - Formatting integers using "sprinti" ; - Changing the line color and thickness in an XY plot ; ;*************************************************************** 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 ;------------------------------------------------------------- ; User specifications ;------------------------------------------------------------- NC = False ; create netCDF? PLOT = True ; sample plots? ymdStrt = 18500101 ; start yyyymmdd ymdLast = 20051231 ; last yrStrt = ymdStrt/10000 yrLast = ymdLast/10000 var = "pr" ; name of file vName = "Daily Anomalies" ; name for plots nhar = 4 ; number of fourier comp ;if (NC) then ; diro= "/home/ahaque/" ; output dir ; filo= var+".day.anomalies."+yrStrt+"-"+yrLast+".nc" ; output file ; end if if (PLOT) then wksType = "x11" wksName = "Daily Anomalies" ; var+"."+yrStrt+"_"+yrLast end if ; Read user specified time and create required yyyyddd ;*********************************************************** a=addfile("/Projects/CMIP5A/Daily/historical/CNRM-CM5/r1i1p1/pr_day_CNRM-CM5_historical_r1i1p1_18500101-20051231.nc","r") pr0= a->pr ;dimensions and size;56978x128x256 ;timexlatxlon time = a->time ; all times on file ymd = cd_calendar(time, -2) ; yyyymmdd iStrt = ind(ymd.eq.ymdStrt) ; index start iLast = ind(ymd.eq.ymdLast) ; index last delete(time) delete(ymd) time = a->time(iStrt:iLast) ; time:units = "hours since" TIME = cd_calendar(time, 0) ; type float year = floattointeger( TIME(:,0) ) month = floattointeger( TIME(:,1) ) day = floattointeger( TIME(:,2) ) ddd = day_of_year(year, month, day) yyyyddd = year*1000 + ddd ; needed for input ;*********************************************************** Convert kg ;mm/day ;*********************************************************** printVarSummary(pr0) printMinMax(pr0,1) pr0 = pr0*86400. pr0@units = "mm/day" ;*********************************************************** ; Compute daily climatology: raw and then 'smoothed' ;*********************************************************** pr0ClmDay = clmDayTLL(pr0, yyyyddd) ; daily climatology at each grid point printVarSummary(pr0ClmDay) ;*********************************************************** ; Compute smoothed daily climatology using 'nhar' harmonics ;*********************************************************** pr0ClmDay_sm = smthClmDayTLL(pr0ClmDay, nhar) printVarSummary(pr0ClmDay_sm) ;*********************************************************** ; Compute daily anomalies using raw and smoothed climatologies ;*********************************************************** pr0Anom = calcDayAnomTLL (pr0, yyyyddd, pr0ClmDay) printVarSummary(pr0Anom) printMinMax(pr0Anom, True) pr0Anom_sm = calcDayAnomTLL (pr0, yyyyddd, pr0ClmDay_sm) pr0Anom_sm@long_name = "Anomalies from Smooth Daily Climatology" printVarSummary(pr0Anom_sm) printMinMax(pr0Anom_sm, True) delete( pr0 ) ; no longer needed ;*********************************************************** ; Create netCDF: convenience use 'simple' method ;*********************************************************** dimx = dimsizes(pr0Anom) ntim = dimx(0) nlat = dimx(1) mlon = dimx(2) if (NC) then lat = a->lat lon = a->lon system("/bin/rm -a "+diro+filo) ; rm any pre-exist file, if any fnc = addfile (diro+filo, "c") filAtt = 0 filAtt@title = vName+": Daily Anomalies: "+yrStrt+"-"+yrLast filAtt@source_file = "original-file.nc" filAtt@creation_date = systemfunc("date") fileattdef( fnc, filAtt ) ; copy file attributes setfileoption(fnc,"DefineMode",True) varNC = vName+"_anom" varNC_sm = vName+"_anom_sm" dimNames = (/"time", "lat", "lon"/) dimSizes = (/ -1 , nlat, mlon/) dimUnlim = (/ True , False, False/) filedimdef(fnc,dimNames,dimSizes,dimUnlim) filevardef(fnc, "time" ,typeof(time),getvardims(time)) filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat)) filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon)) filevardef(fnc, varNC ,typeof(pr0Anom) ,getvardims(pr0Anom)) filevardef(fnc, varNC_sm,typeof(pr0Anom) ,getvardims(pr0Anom)) filevarattdef(fnc,"time" ,time) ; copy time attributes filevarattdef(fnc,"lat" ,lat) ; copy lat attributes filevarattdef(fnc,"lon" ,lon) ; copy lon attribut filevarattdef(fnc,varNC ,pr0Anom) filevarattdef(fnc,varNC_sm,pr0Anom_sm) fnc->time = (/time/) fnc->lat = (/lat/) fnc->lon = (/lon/) fnc->$varNC$ = (/pr0Anom /) fnc->$varNC_sm$ = (/pr0Anom_sm/) end if ;************************************************ ;Plotting parameters and resources ;************************************************ ;Plotting parameters and resources ;************************************************ ; Plot anomalies for an arbitrarily selected near equatorial location ; Time: Oct 1, 1996 to April 1,1997 [arbitrary selection] ;========== LATX = 26.9 ; lat for south asia LONX = 72.2 ; lon for south asia yyyymmdd = cd_calendar(time, -2) pr0Anom@long_name = "Anomalies from Raw" ; short labels for plot pr0Anom_sm@long_name = "Anomalies from Smooth" ntBegin = ind(yyyymmdd.eq.18500101) ntEnd = ind(yyyymmdd.eq.20051231) monthLabels = (/1,4,7,10/) monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \ ,"Jul","Aug","Sep","Oct","Nov","Dec" /) pr0Val = new(ntim, typeof(pr0Anom&time) , "No_FillValue") ; bigger than pr0Lab = new(ntim, "string", "No_FillValue") ; needed pr0Valm = new(ntim, typeof(pr0Anom&time) , "No_FillValue") ; bigger than ntm = -1 cr = inttochar(10) ; carriage return do nt=ntBegin,ntEnd if (day(nt).eq.1) then ntm = ntm + 1 pr0Val(ntm) = pr0Anom&time(nt) pr0Lab(ntm) = monNam(month(nt)-1) if (month(nt).eq.1) then pr0Lab(ntm) = pr0Lab(ntm) + cr +sprinti("%0.4i", year(nt)) end if end if end do wks = gsn_open_wks ("x11"," Daily Anomalies") rxy = True rxy@gsnDraw = False rxy@gsnFrame = False rxy@gsnYRefLine = 0.0 ; create a reference line rxy@gsnAboveYRefLineColor = "red" ; above ref line fill red rxy@gsnBelowYRefLineColor = "blue" ; below ref line fill blue rxy@xyLineThicknessF = 2.0 rxy@vpHeightF = 0.4 ; resize rxy@vpWidthF = 0.8 rxy@tmXBMode = "Explicit" rxy@tmXBValues = pr0Val(0:ntm) rxy@tmXBLabels = pr0Lab(0:ntm) plot(0) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \ ,pr0Anom(ntBegin:ntEnd,{26.9},{72.2}),rxy) plot(1) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \ ,pr0Anom_sm(ntBegin:ntEnd,{26.9},{72.2}),rxy) resP@txString = vName+": Anomalies: (26.9,72.2E)" gsn_panel(wks,plot(0:1),(/2,1/),resP) end