##COLLATE A1B MODEL RUNS FOR RAHMSTORF ANALYSIS

###FUNCTIONS 
	source("http://www.climateaudit.org/scripts/models/functions.collation.knmi.txt")	
		#function to scrape from KNMI
	
## LOGON
	email=Email=   ##register and insert your email here
	logon=download_html( paste("http://climexp.knmi.nl/start.cgi?",Email,sep=""))
  	scenario="sresa1b"
	field="tas";
	lat0=c(-90,90)
	Info=knmi.info[knmi.info$scenario==scenario,] 
	row.names(Info)=1:nrow(Info)
	Info #gives A1B models at KNMI 

	#              model             alias scenario Runs
	#1            BCC CM1           bcc_cm1    20c3m    4
	#2        BCCR BCM2.0       bccr_bcm2_0    20c3m    1
	#4      CGCM3.1 (T47)     cccma_cgcm3_1    20c3m    5
	

	tempsc=knmi.info$scenario==scenario;sum(tempsc);M=sum(tempsc)
	   #select models with A1B runs

#scrape 24 models from KNMI (July 2009)
	Model=rep(list(NA),M);names(Model)=knmi.info$alias[tempsc] #placeholder 

	for (i in 1:sum(tempsc) ) {
	  model=knmi.info[tempsc,][i,"alias"]
   	  Model[[i]]= read.knmi.models(field="tas",model ,scenario, lat=lat0,version="raw") 
	       #scrapes all runs for this model into time series
	 }

#clean up NA values  
	for( i in 1:length(Model) ) { temp= Model[[i]]< -999|Model[[i]] > 1E6 
		Model[[i]][temp]=NA
		}

#inspection of HadCM3 showed weird March values. One run excluded
	sapply(Model,function(A) range(A,na.rm=T) )
	  #this shows weird 0 values for HadCM3 (only one run)
	Model=Model[c(1:22,24)] #exclude HAD problem

#require all models to have 1961-1990 values for reference period
	(test=sapply(Model, function(A) tsp(A)[1]) )
      		# bccr_bcm2_0     cccma_cgcm3_1 cccma_cgcm3_1_t63          cnrm_cm3       csiro_mk3_0       csiro_mk3_5        gfdl_cm2_0 
      		#       1850              1850              1850              1860              1871              1871              1861 
      		# gfdl_cm2_1          giss_aom    giss_model_e_h    giss_model_e_r   iap_fgoals1_0_g       ingv_echam4          inmcm3_0 
      		#       1861              2001              1880              1880              2000              1870              1871 
  
	temp=test> 1960 ;sum(!temp) #21
	Model=Model[!temp];length(Model) #21	

#window each set of model runs to 1900-2099.99
	Model= sapply(Model, function(A) window(A,start=1900,end=2099.99) )

#make anomaly with rference 1961-1990 
	Anom=sapply(Model, anomaly )

#convert R-list into time series with 51 columns 
	name0=unlist( sapply(Anom,function(A) dimnames(A)[[2]]) );name0
 		#     bccr_bcm2_0        cccma_cgcm3_11        cccma_cgcm3_12 ...
      		#bccr_bcm2_0_1"     "cccma_cgcm3_1_1"     "cccma_cgcm3_1_2" ...

	working=NULL
	for (i in 1:length(Anom)) working=ts.union(working,Anom[[i]])
	for(i in 1:ncol(working)) working[,i]=anomaly(working[,i])
	dim(working) #2400 51

#exclude runs with no 1961-1990 values
	temp= time(working)>=1961 & time(working)<1991
	m0=apply(working[temp,],2,mean);m0
	temp=is.na(m0);sum(temp) #2 PCM models without 1961-1990 values

	working=working[,!temp];dim(working)

	year=c(t(  array( rep(1900:2099, 12),dim=c(200,12) ) ))
	month= rep(1:12,200)
	working=cbind(year,month,working)
	dimnames(working)[[2]]=c("year","month",name0[!temp])

	zz <- gzfile("d:/climate/data/rahmstorf/a1b_models.gz", "w")  # compressed file
	write.table(working,zz,sep="\t",row.names=FALSE)
	close(zz)
	list.files("d:/climate/data/rahmstorf")

	url= "d:/climate/data/rahmstorf/a1b_models.gz"
	handle <- gzfile(url) 
	working= read.table(handle,sep="\t",header=TRUE)
	close(handle)
	list.files("d:/climate/data/rahmstorf")

		dim(working) #2400 53

	#this is uploaded to www.climateaudit.org/data/rahmstorf/a1b_models.gz 











readLines(zz <- gzfile("ex.gz"))
close(zz)
unlink("ex.gz")



 handle=gzfile(