### SANTER FUNCTIONS ### ### df.d1star=function(sd.model,M,se.obs,neff) ### Santer equation 13 made into a function to calculate df for t test ### ### santerf - calculates trend significance using Santer AR1 method ### ### function(x,end0=2050,start0=1979,M=19, ### trend.target=trend.model.T2LT, sd.target=sd.model.T2LT ) ### ### ### ### ### ### #### ### ### ### ### ### ### ### ####### df.d1star df.d1star=function(sd.model,M,se.obs,neff) { #df calculation from Santer equation 13 made into function C= (sd.model/sqrt(M-1))^2 + se.obs^2;C # 0.0195 D=(sd.model/sqrt(M-1))^2 /(M-1) + se.obs^2/(neff-2) ;D # 0.00155 df.d1star= C/D # 12.6 df.d1star} #### santerf santerf=function(x,end0=2050,start0=1979,M=19, trend.model=trend.model.T2LT, sd.model=sd.model.T2LT ) { #f=rep(NA,9); names(f)=c("trend","se","sd","r1","neff","adj_df", "dstar","d1star","d1zero","pt","ptzero") #Douglass stat not used other than in benchmarking f=rep(NA,10); names(f)=c("trend","se","sd","r1","neff","adj_df", "d1star","d1zero","pt","ptzero") temp= (time(x)< (end0+1))&(time(x)>=start0) x=ts(x[temp],start=c(start0,1),freq=12) year=c(time(x))/10; N=length(x); year=year-mean(year) f["sd"]=sd(x) fm= lm (x~year) trend.obs=f["trend"]= fm$coef[2] r=f["r1"]= arima (fm$residuals,order=c(1,0,0))$coef[1]; # neff=f["neff"]= N * (1-r)/(1+r) ; se.obs=f["se"] = sqrt((N-2)/(neff-2))* summary(fm)$coef[2,"Std. Error"] f["d1star"]= (trend.model-trend.obs)/ sqrt( (sd.model/sqrt(M-1))^2 + se.obs^2) ##f["dstar"]= (trend.model-trend.obs)/ (sd.model/sqrt(M-1)) f["adj_df"]= df.d1star(sd.model,M,se.obs,neff) f["d1zero"]=f["trend"]/f["se"] f["pt"]=pt(f["d1star"],f["adj_df"]) f["ptzero"]=pt(f["d1zero"],f["adj_df"]) santerf=f santerf } ###LOAD DATA #1. Santer 2008 Table 1 make.table1=function(loc="http://www.climateaudit.org/data/models/santer_2008_table1.dat") { santer=read.table(loc,skip=1) names(santer)=c("item","layer","trend","se","sd","r1","neff") row.names(santer)=paste(santer[,1],santer[,2],sep="_") santer=santer[,3:ncol(santer)] #insert Table III info santer$d1star=NA;santer$dstar=NA santer["UAH_T2LT",c("d1star","dstar")]= c(1.11,7.16) #from Table III for UAH santer["RSS_T2LT",c("d1star","dstar")]= c(0.37,2.25) #from Table III for UAH santer["UAH_T2",c("d1star","dstar")]= c(1.19,6.78) #from Table III for UAH santer["RSS_T2",c("d1star","dstar")]= c(0.44,2.48) #from Table III for UAH make.table1= santer make.table1 } # trend se sd r1 neff d1star dstar #HadCRUT3v_TL+O 0.119 0.117 0.197 0.934 8.6 NA NA #Multi-model_mean_TL+O 0.146 0.214 0.274 0.915 11.7 NA NA #... make.table4=function() { make.table4=read.csv("d:/climate/data/models/santer/santer_2008_table4.csv",sep="\t",header=TRUE) #this also has Table VI info make.table4 } table1=santer=make.table1("d:/climate/data/models/santer/santer_2008_table1.dat") table4=make.table4() tropo=c("T0","T2LT","T2"); firma=c("All","Ocean","Land") p=length(tropo);q=length(firma) make.model.info=function() { model.info=data.frame(tropo=gl(p+1,q,labels=c(tropo,"Lapse")),firma=rep(firma,p+1),trend=rep(NA,(p+1)*q),sd=rep(NA,(p+1)*q)) model.info[1,3:4]=table1$trend[2:3] #T0 All from S08 T1 model.info[2,3:4]=table1$trend[7:8] #T0 Ocean from S08 T1 model.info[4,3:4]=table1$trend[11:12] #T2LT All from S08 T1 model.info[7,3:4]=table1$trend[15:16] #T2 All from S08 T1 model.info[10,3:4]=table4$Trend[3:4] #T2LT All Lapse from S08 T4 model.info[11,3:4]=table4$Trend[11:12] #T2LT Ocean Lapse from S08 T4 model.info$trend[5]=model.info$trend[2]-model.info$trend[11] #T2LT Ocean= T0 Ocean minus Lapse Ocean-T2LT k= 0.7822222 #proportion TRP ocean from hadISST mask model.info$trend[3]=(model.info$trend[1]- k*model.info$trend[2])/(1-k) #T0 Land: by difference from T0 All and T0 Ocean model.info$trend[6]=(model.info$trend[4]- k*model.info$trend[5])/(1-k) #T2LT by difference from T2LT All and T2LT Ocean: no difference model.info$trend[12]=(model.info$trend[3]- model.info$trend[6]) #T2LT Land Lapse: T0 Land minus T2LT Land model.info$sd[5:6]=model.info$sd[4] ###interim estimate since not reported by Santer model.info$sd[3]=.07 #interim estimate higher than T0 Ocean model.info$sd[12]=.04 #estimate interim higher than Lapse Ocean return(model.info) } model.info=make.model.info() tropo firma trend sd 1 T0 All 0.14600000 0.066 2 T0 Ocean 0.13000000 0.062 3 T0 Land 0.20346938 0.070 4 T2LT All 0.21500000 0.092 5 T2LT Ocean 0.21500000 0.092 6 T2LT Land 0.21500000 0.092 7 T2 All 0.19900000 0.098 8 T2 Ocean NA NA 9 T2 Land NA NA 10 Lapse All -0.06900000 0.032 11 Lapse Ocean -0.08500000 0.038 12 Lapse Land -0.01153062 0.040 ### #msy[,"Trpcs_ocean"] source("http://www.climateaudit.org/scripts/spaghetti/msu.glb.txt") source("http://www.climateaudit.org/scripts/spaghetti/msu.t2.txt") rss=list(tropo=c("TLT","TMT"),firma=c("Land_and_Ocean","Ocean","Land")) get.rss=function(tropo,firma) { ##tropo="TLT";firma="Ocean" loc=paste("http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_",tropo,"_Anomalies_",firma,"_v03_2.txt",sep="") get.rss=read.table(loc, skip=3) name0=c("year","month","82.82","20.20","20.82N","20.82S","60.80N","60.80S") if(!(firma=="Land_and_Ocean")) names(get.rss)=name0 else names(get.rss)=c(name0,"US","NH","SH") get.rss=ts(get.rss[,"20.20"],start=1979,freq=12) get.rss } get.uah=function(tropo,firma) { if(firma=="All") { if(tropo=="T2LT") get.uah=msu[,"Trpcs"] if(tropo=="T2") get.uah=msu.t2[,"Trpcs"] } else { if(tropo=="T2LT") get.uah=msu[,paste("Trpcs",casefold(firma),sep="_")] if(tropo=="T2") get.uah=msu.t2[,paste("Trpcs",casefold(firma),sep="_")] } get.uah } do.santer=function( X, target, method="lapse") { if ( method=="level") { Obs=rep(list(NA),3);names(Obs)=c(1999,2007,2009); q=ncol(X) for (k in 1:3) Obs[[k]]=data.frame( id=dimnames(X)[[2]],target=rep( target["trend"],q),trend=NA,se=NA,sd=NA,r1=NA,neff=NA,adj_df=NA,d1star=NA,d1zero=NA,pt=NA,ptzero=NA) for (j in 1:q) { temp1=!is.na(X[,paste(Obs[[1]]$id[j])] ); end0=min( max(time(X)[temp1]) ) x=window( X[,paste(Obs[[1]]$id[j])],start=1979,end=end0 ) Obs[[1]][j,3:12]=round(f(x,end0=1999,trend.model=target["trend"] ,sd.model=target["sd"]),3) Obs[[2]][j,3:12]=round(f(x,end0=2007,trend.model=target["trend"] ,sd.model=target["sd"]),3) Obs[[3]][j,3:12]=round(f(x,end0=2009,trend.model=target["trend"] ,sd.model=target["sd"]),3) } for(k in 1:3) row.names(Obs[[k]])=paste("#",row.names(Obs[[k]])) return(Obs) } else { Surf=X$surf; Sat=X$sat; Lapse=rep(list(NA),3);names(Lapse)=c(1999,2007,2009) p=ncol(Surf); q=ncol(Sat) for (k in 1:3) { Lapse[[k]]=data.frame( surf=gl(p,q,labels=dimnames(Surf)[[2]]), sat=rep(dimnames(Sat)[[2]],p),target=rep( target["trend"],p*q), trend=NA,se=NA,sd=NA,r1=NA,neff=NA,adj_df=NA,d1star=NA,d1zero=NA,pt=NA,ptzero=NA) } #santerf= function(x,end0=2050,start0=1979,M=19, #trend.model=trend.model.T2LT, sd.model=sd.model.T2LT ) { for (j in 1:(p*q)) { temp1=!is.na(Surf[,paste(Lapse[[1]]$surf[j])] ); temp2=!is.na(Sat[,paste(Lapse[[1]]$sat[j])] ) ; end0=min( max(time(Surf)[temp1]), max(time(Sat)[temp2]) ) x=window( Surf[,paste(Lapse[[1]]$surf[j])],start=1979,end=end0 ) y=window( Sat[,paste(Lapse[[1]]$sat[j])],start=1979,end=end0 ) lapse=x-y Lapse[[1]][j,4:13]=round(f(lapse,end0=1999,trend.model=target[1],sd.model=target[2]),3) Lapse[[2]][j,4:13]=round(f(lapse,end0=2007,trend.model=target[1],sd.model=target[2]),3) Lapse[[3]][j,4:13]=round(f(lapse,end0=2009,trend.model=target[1],sd.model=target[2]),3) } for(k in 1:3) row.names(Lapse[[k]])=paste("#",row.names(Lapse[[k]])) return(Lapse) } } ##MODEL COLLATION FUNCTIONS #load("d:/climate/data/models/santer/tam.tab") info.pcmdi=read.csv("d:/climate/data/models/santer/info.pcmdi.csv",sep="\t",header=TRUE) #collated in collation.pcmdi.txt info.model=read.csv("d:/climate/data/models/santer/info.model.csv",sep="\t",header=TRUE) #collated in collation.pcmdi.txt #id model runs alias #1 1 CCCMA3.1 1 CGCM3.1 (T63) info.model=info.model[,1:4] volc=c(0,1,0,0,0,1,1,0,1,1,0,1,0,0,0,1,1,1,1) #compiled from Santer PNAS Figure 4d info.model$volc=factor(volc); levels(info.model$volc)=c("NV","V") ##COLLATE T2 and T2LT Time Series make.santermodel=function(Tam=tam,Info.pcmdi=info.pcmdi) { #info.pcmdi=read.csv("http://www.climateaudit.org/data/models/santer/info.pcmdi.csv",sep="\t",header=TRUE) temp_model=!is.na(Info.pcmdi$id);sum(temp_model) #48 test=sapply(Tam[["T2"]],function(A) A[,"TROP"]) T2=NULL for(i in 1:length(test)) T2=ts.union(T2,test[[i]]) T2=T2[,temp_model] dimnames(T2)[[2]]=Info.pcmdi$id[temp_model] temp= T2< -900 T2[temp]=NA #dim(T2) #[1] 2412 48 #tsp(T2) # 1850.000 2050.917 12.000 test=sapply(Tam[["T2LT"]],function(A) A[,"TROP"]) T2LT=NULL for(i in 1:length(test)) T2LT=ts.union(T2LT,test[[i]]) T2LT=T2LT[,temp_model] dimnames(T2LT)[[2]]=Info.pcmdi$id[temp_model] temp= T2LT< -900 T2LT[temp]=NA #dim(T2LT) #[1] 2412 49 #tsp(T2LT) # make.santermodel=list(T2LT=T2LT,T2=T2) make.santermodel } #SanterModel=make.santermodel()