##SCRIPT FOR CALCULATIONS IN SANTER COMMENT


########################
##LOAD FUNCTIONS AND DATA
###########################

library(nlme)
source("http://www.climateaudit.info/scripts/utilities.txt")
source("http://www.climateaudit.info/scripts/models/santer/functions.santer.txt") ; 
  # returns: table1; table4; santer.model.info; plus some functions 

################################
## Load Observations
#############################
  # these have been previously collated as at March 10, 2010 by v3/collation.observations.txt

  ## sat Observations
	download.file("http://www.climateaudit.info/data/models/santer/2010/Sat_201004.tab","temp.dat",mode="wb"); load("temp.dat") #Sat
	#load("d:/climate/data/models/santer/Sat_201004.tab")
	   #collated in santer/v3/collation.observations.txt
 	 Sat[["T2LT"]][["All"]][1,]
	 #    rss    uah 
	 # -0.249 -0.090 

  ## SURF Observations
	download.file("http://www.climateaudit.info/data/models/santer/2010/Surf_201004.tab","temp.dat",mode="wb"); load("temp.dat") #Sat
	#load("d:/climate/data/models/santer/Surf_201004.tab")
	   #Surf: All, OCean, Land 
	   #collated in santer/v3/collation.observations.txt

  ## Radiosonde Observations
	download.file("http://www.climateaudit.info/data/models/santer/Radio_201002.tab","temp.dat",mode="wb"); load("temp.dat") #Sat
	#load("d:/climate/data/models/santer/Radio_201005.tab")


#################################
## Load Runs and Run Information
####################################

	download.file("http://www.climateaudit.info/data/models/chad/info.chad.tab","temp.dat",mode="wb"); load("temp.dat") #Sat
	# load("d:/climate/data/models/chad/info.chad.tab")
	names(info.chad) #  "models" "runs"  
	  info.chad$runs[1:2,]
		#                       id              model run              santer
		#1       bccr_bcm_2.0_run_1       bccr_bcm_2.0   1                <NA>
		#2 cccma_cgcm_3.1_T47_run_1 cccma_cgcm_3.1_T47   1 CCCMA3.1_20c3m_run1
      ## don't use miub_echo_g which lacks atmospheric
	temp = info.chad$runs$model=="miub_echo_g" ; sum(!temp)  #82
      #restrict to A1B 
	temp1=info.chad$runs$end>2004
	Info=info.chad$runs[!temp&temp1,1:3] 
	dim(Info) #57 3
       	Info$model=factor(Info$model)
 	length(levels(Info$model)) 	# 24
	Info.chad=Info

#load runs
	download.file("http://www.climateaudit.info/data/models/chad/Runs.anomaly.tab","temp.dat",mode="wb"); load("temp.dat") #Sat
	#load("d:/climate/data/models/chad/Runs.anomaly.tab") 
	   #v4/collation.averages.txt using v4/collation.models.txt
	Runs=Runs.anomaly
	names(Runs) # 	"land"  "ocean" "loti" 
	   sapply(Runs,length)
 	   	#land ocean  loti 
    	   	#3     3     3 
	   sapply(Runs$loti,tsp)
	   	#        T0      T2LT      T2
		#[1,] 1850.000 1850.000 1850.000
		#[2,] 2099.917 2099.917 2099.917
	   sapply(Runs$loti,dim)
 		#     TAS  TLT  TMT
		#[1,] 3000 3000 3000
		#[2,]   57   57   57

#################
#calculate target info for models ended in 2009
#######################
    #this is a utility function that collects trend and sd information from models for use in 
    #santer-style comparisons
	Data.A1B= rep(list(NA),3) ;names(Data.A1B)= c("land","ocean","loti")
	for (i in 1:3) Data.A1B[[i]]= collate.model.information(Runs=Runs.anomaly[[i]],end0=2009.99)
	Data.A1B$loti$target.info 

		#         trend    sd count
		#T0       0.193 0.060    24
		#T2LT     0.272 0.082    24
		#T2       0.262 0.087    24
		#T0_T2LT -0.079 0.029    24
		#T0_T2   -0.069 0.034    24


####################################
## CALCULATE TRENDS WITH SANTER-STYLE AR1 SIGNIFICANCE
####################################

  #uses utility function to collect trend information for variety of cases
  #these are santer comparisons

	Trend.A1B=make.trend.information(Data.A1B)

	Trend.A1B$Lapse_T2LT$All

	  #     obs target d1star    pt santer
	  #rss -0.033 -0.079 -2.819 0.003   0.67
	  #uah  0.048 -0.079 -7.395 0.000   3.50

	Trend.A1B$Lapse_T2$All

	   #       obs target d1star pt santer
	   #rss 0.005 -0.069 -4.212  0     NA
	   #uah 0.084 -0.069 -8.518  0     NA


	Trend.A1B$T2LT$All

	  #       obs target d1star    pt santer
	  #rss 0.159  0.272  1.690 0.948   0.37
	  #uah 0.075  0.272  2.862 0.996   1.11

	Trend.A1B$T2$All

	  #          obs target d1star    pt santer
	  #rss 0.121  0.262  2.196 0.981   0.44
	  #uah 0.040  0.262  3.449 0.999   1.19

	
	Trend.A1B$T0$All
 
     	  #       obs target d1star    pt santer
	  #had  0.123  0.193  1.179 0.870     NA
	  #giss 0.115  0.193  1.333 0.898     NA
	  #noaa 0.124  0.193  1.327 0.897     NA