#Load in the climate RCPs
download.file("http://www.stat.washington.edu/peter/593/ForPeter/glbtas.cmip5.rcps.1850-2300.r1i1p1.nc","glbtas.cmip5.rcps.1850-2300.r1i1p1.nc")

cmip5=open.ncdf("glbtas.cmip5.rcps.1850-2300.r1i1p1.nc")
cmip5.data=get.var.ncdf(cmip5,"glbtas") #global mean temperatures

#turn 9e36 values into NAs
for(i in 1:451) {for(j in 1:45) {for(k in 1:4) cmip5.data[i,j,k]=ifelse(cmip5.data[i,j,k]>9e36,NA,cmip5.data[i,j,k])}}

cmip5.data=cmip5.data-273.15 #Convert from Kelvin to °C

for(i in 1:45) {for (k in 1:4) cmip5.data[,i,k] = cmip5.data[,i,k]-mean(cmip5.data[121:150,i,k],na.rm=T)} #anomalies wrt 1970-1999

#pick only models with complete data 2000-2100, should return 18 models
models=NULL
for(i in 1:45) if(sum(sum(is.na(cmip5.data[151:251,i,])))==0)models=c(models,i)