#RAOBCORE VERSION 1.4 #see commentary in raobcore.txt source("http://www.climateaudit.info/scripts/utilities/getx.var.ncdf.txt") get.raobcore.trp=function(loc="ftp://srvx6.img.univie.ac.at/pub/raobcore14_gridded_2009.nc",method="trp") { download.file(loc,"zonal.nc",mode="wb") v<-open.ncdf("zonal.nc") v$var[[1]]$dim[2][[1]]$vals #anomalies # 85 75 65 55 45 35 25 15 5 -5 -15 -25 -35 -45 -55 -65 -75 -85 instr <- getx.var.ncdf( v, v$var[[1]]) # 1850 2006 dim(instr) # 36 18 12 612 temp= abs(v$var[[1]]$dim[2][[1]]$vals)<20;sum(temp) #4 instr=instr[,temp,,];dim(instr) # 36 4 12 612 instr=aperm(instr,c(4,1,2,3)); dim(instr) # 588 36 4 12 dim0=dim(instr) instr=array(instr,dim=c(dim0[1],144,12) ); dim(instr) # 612 144 12 raobcore= ts( array (NA,dim=c(dim0[1],12) ),start= 1958,freq=12); tsp(raobcore) # 1958.000 2008.917 12.000 for(i in 1:12) raobcore[,i]= apply(instr[,,i],1,mean,na.rm=T) dimnames(raobcore)[[2]]=v$var[[1]]$dim[3][[1]]$vals # 850 700 500 400 300 250 200 150 100 70 50 30 if(method=="T2LT") { christy.weights= c(0.1426, 0.0706, 0.1498, 0.1894, 0.1548, 0.1382, 0.1035, 0.0477, 0.0167, 0.0067, -0.0016, -0.0179) ##CHRISTY WEIGHTS email Apr 30 2008 names(christy.weights)= c(1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100) raobcore_t2lt= apply(raobcore[,c("850","700")],1, function(x) weighted.mean(x,c(.33,.67)) ) raobcore_t2lt=ts (raobcore_t2lt,start=1958,freq=12) raobcore_t2lt=round(raobcore_t2lt,3) raobcore=raobcore_t2lt} return(raobcore) } #raobcore=get.raobcore.trp() #. Use Christy Weights to make T2LT weighted average of trends #effective T2LT altitude is 2.5 km, about 750 hPa #christy.weights= c(0.1426, 0.0706, 0.1498, 0.1894, 0.1548, 0.1382, 0.1035, # 0.0477, 0.0167, 0.0067, -0.0016, -0.0179) ##CHRISTY WEIGHTS email Apr 30 2008 #names(christy.weights)= c(1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100) #T2 - mean altitude 6.1 km) 466 hPa #http://www.sensorsone.co.uk/altitude-pressure-units-conversion.html #effective T2LT altitude is 2.5 km, about 750 hPa (Christy) #TLS is centered around 75 hPa (~20 km or ~14-50 km) and TTS (channel 3) around 220 hPa (~10 km or 1-35 km). #. Use Christy Weights to make T2LT weighted average of trends #effective T2LT altitude is 2.5 km, about 750 hPa #christy.weights= c(0.1426, 0.0706, 0.1498, 0.1894, 0.1548, 0.1382, 0.1035, # 0.0477, 0.0167, 0.0067, -0.0016, -0.0179) ##CHRISTY WEIGHTS email Apr 30 2008 # names(christy.weights)= c(1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100) # fred=apply(window(raobcore,start=1979),2,trend) # 850 700 500 400 300 250 200 150 100 70 50 30 # 0.195 0.042 0.032 0.132 0.280 0.295 0.202 0.009 -0.151 -0.477 -0.567 -0.613 # plot(fred, as.numeric(names(fred)),log="y",ylim=c(1000,100),yaxs="i",xlim=c(-.6,.5),col=4,type="b") # points( trend(Sat$T2LT$All[,"uah"]) ,750,pch=19) # points( trend(Sat$T2$All[,"uah"]) ,466,pch=19) # points( trend(Sat$T2LT$All[,"rss"]) ,750,pch=19,col=2) ## points( trend(window(Surf$All[,"had"],start=1979)) ,999,pch=16,col=3,cex=2) # fred1=apply( window(hadat,start=1979),2,trend) # id=names(fred1);n=nchar(id); names(fred1)=substr(names(fred1),1,n-3) # points(fred1,as.numeric(names(fred1)),col=2,type="b",pch=18) # ratpac=ratpac[,2:14] # fred1=apply( window(ratpac,start=1979),2,trend) # names(fred1)[1]=1000; # points(fred1/12,as.numeric(names(fred1)),col=3,type="b",pch=18)